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I.  INTRODUCTION 


1.1  DETECTION  THEORY 

The  problem  of  the  detection  of  a  seismic  signal,  sometimes 
present,  sometimes  absent,  In  the  presence  of  noise  may  be  conveniently 
discussed  using  the  theory  z>  statistical  detection.  {See  Helstrom, 

1968,  as  a  good  general  reference.)  Based  on  a  set  of  measurements,  we 
are  trying  to  choose  between  two  hypotheses: 

Hq  that  the  measurements  consist  only  of  noise 

that  the  measurements  consist  of  signal  plus  noise. 

For  discrete  data,  we  can  conveniently  use  matrix  notation  and  represent 
the  set  of  measurements  as  a  vector  x.  with  cu.iipoenents  that  are 
the  individual  samples  of  the  seismometer's  output  (the  detector's  in¬ 
put). 

The  number  of  components,  N,  in  this  vector  is  the  number  of  samples 
of  the  input  that  we  consider  (process)  at  one  time.  We  can  measure  some¬ 
thing  about  the  characteristics  or  statistics  of  x  over  a  given  time 
interval,  and  a  signal  may  be  detected  if,  over  this  interval,  these  in¬ 
put  statistics  are  significantly  different  from  what  we  would  expect 
from  noise  alone.  Now,  obviously,  we  will  maximize  this  statistical 
difference  if  we  match  this  time  interval  to  the  finite  duration  of  the 
expected  seismic  signal.  This  fundamental  time  duration,  TD,  is  typ¬ 
ically  on  the  order  of  a  second,  and  a  quantity  calculated  over  this 
interval  is  often  referred  to  as  the  "short  term  average"  (STA)  of  the 
input  signal.  Thus  the  number  of  components  in  the  vector  ^  will  be 


where  At  is  the  sampling  interval. 
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The  detector  must  decide  between  the  hypothesis  Hq  and  Ay 
The  correct  strategy  to  use  depends  upon  the  cost  or  the  risk  of 
making  the  wrong  decision.  In  seismic  applications  (as  opposed  to  game 
theory  for  example)  the  cost  or  risk  of  making  the  wrong  decision  is 
Impossible  (or  at  least  difficult)  to  quantify.  Instead,  what  is 
commonly  done  Is  to  use  the  Neyman-Pearson  strategy  which  maximizes 
the  probability  of  detection  with  a  specified  probability  of  being  wrong  - 
a  false  alarm.  Conceptually,  we  can  think  of  the  process  of  hypothesis 
testing  by  using  the  probability  density  functions  PqOO  and  P-jU) 
of  Hq  and 


Figure  1.1.  Probability  density  function  of  the  Hq  and 
H1  hypotheses. 

The  chance  of  error,  Qq,  caused  by  choosing  when  Hq  is  true, 
that  is  "detecting"  a  signal  when  one  is  not  there,  is  called  a  false 
alarm  probability.  It  is  given  by 

oo 

Qq  -/  Po^)dI  (1.1) 

*0 
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Observe  that: 


1.  The  probability  of  false  alarm,  Qq,  and  the  false  alarm 
rate,  FAR  (number  of  false  alarms  per  unit  time),  are 
related  by  a  time  Interval  Tq 


where  Tq  Is  the  time  Interval  (non  overlapping)  over  which 
the  "statistic"  Is  calculated. 

2.  Increasing  Tq  decreases  FAR  for  a  fi  <ed  Qq  but  also  de¬ 
creases  the  time  resolution  of  the  detector. 

3.  Increasing  TD  also  provides  for  more  frequency  resolution 
for  those  detectors  working  in  the  frequency  domain. 

4.  In  the  case  where  there  is  a  signal  present,  increasing  Tq 
will  possibly  decrease  the  signal-to-noise  ratio  (SNR)  by 
including  more  noise  after  the  signal  has  died  away. 

The  error,  Qj,  of  choosing  HQ  when  Hj  is  true,  that  is,  not 
detecting  a  signal  when  one  is  present,  is  given  by 


(1.2) 


and  the  detection  probability  Qq  is 

Qo  ■  1  -  Qx  •  (1.3) 

The  Neyman-Pearson  strategy  of  decision  theory  leads  us  to  search  for 
the  minimum  Qj  for  a  specified  Qq.  Ideally,  for  each  detector  we  would 
like  to  be  able  to  construct  a  diagram  like  Figure  1.2  which  in  detection 
theory  is  called  the  "receiver  operating  characteristic"  (ROC). 


3 


SYSTEMS.  SCIENCE  ANO  SOETWAmC 


i  in**-  -**  ’ 


0  i 


% 

FALSE  ALARM  PROBABILITY 

Figure  1.2.  The  receiver  operating  characteristic  (ROC). 


The  slope  of  the  ROC  curve  is 

dQ0  P1 (*) 

W0  =  PjxJ 


which  is  called  the  "likelihood  rat<o".  This  is  easy  to  see  as 


-  H) 

Q0  »  1  -  Q1  ■  f  P1(x)dx  -  J  P-j (^)dx.  *  f  P^xjdx 

_  u 


ao 


i.(i)Po(i>dx 


(1.5) 


Thus 


dQ 


dQ, 


■af "  a(*>p0(*>  snd  -ar  ■  p0(i) 


SO 


dQn 

dC  ’  A!i) 


(1.7) 


For  example,  in  the  case  of  a  Gaussian  signal  in  Gaussian 


noise, 


sFoexp(-^M 

(1.8) 

(1.9) 

where  RQ  and  are  the  covariance  matrices  of  noise  and  signal  plus 

noise  and  x*  is  the  transpose  of  x.  We  see  that  the  likelihood  function 

A(x)  -  k  exp||  x^Rq1  -  R^)x) 

(1.10) 

is  a  monotonic  function  of  the  quadratic  form 


x*  A  x 

where 


(1.11) 


(1.12) 


As  He! strom  (1968,  p.  94)  shows,  the  choice  between  hypotheses  HQ 
and  can  equally  be  based  on  any  monotonic  function  of  the  likelihood 
function  A(x).  The  quantity  that  summarizes  and  replaces  that 
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data  x  themselves  Is  tenned  the  detection  "statistic".  Clearly,  then, 
we  could  equally  well  use  the  quadratic  form 

xlAx 

as  our  detection  statistic. 

Me  may  generate  this  statistic  by  the  action  of  a  linear  filter 
H  on  the  signal  followed  by  a  device  that  forms  the  product 

(H  x)t  •  (H  x) 
since  this  is 

X*  M*  H  X  -  x1  A  X  (1.13) 

It  may  be  shown  that  if  A  Is  a  positive  definite  synmetric 

matrix, 

(i.i4) 

where  D  is  the  diagonal  matrix  of  the  square  roots  of  the  eigenvalues 
of  A  while  Q  is  the  matrix  of  normalized  eigenvalues  of  A. 
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1.2  THE  FRIEBERGER  DETECTOR 


In  an  elegant  work  by  Frieberger  (1963)  it  is  shown  that  the 
quadratic  form  given  by  Eq.  (1.11)  can  be  approximated  by  a  constant  times 
the  Integral 


F$(w) 

Fn(w)[Fs(w)  + 


(1.15) 


In  the  case  where  the  signal  and  noise  both  have  a  Gaussian  probability 
density  function.  In  this  equation,  Fx,  F^  and  F$  are  the  power  spectral 
densities  of  the  input,  the  noise  alone,  and  the  signal  alone,  respectively. 
This  detector  can  thus  be  implemented  as  illustrated  in  Figure  1.3. 


x(t) 


>C 

'<C 


H, 


Figure  1.3.  Implementation  of  the  Frieberger  detector. 


where  the  modulus  of  the  frequency  response  of  the  adaptive  filter  is  given 
by: 


V^)[FSU)  "fnU)]  u  j 

and  the  power  meter  is  a  quadratic  integrator.  Because  of  the  squaring  op¬ 
eration,  the  filter  phase  response  is  irrelevant.  Note  that  for  a  large  SNR 
(i.e.,  F$(oj)  »  Fn(w)),  this  response  is  approximately^/l/F^U) .  It  is 
not  the  Weiner  filter  which  is^F$(w)/(Fs(u))  +  FN(w)),  whose  task  Is  not 
to  detect  but  rather  to  recover  the  signal. 

What  has  become  to  be  referred  to  as  the  "conventional  power 
detector"  simply  performs  the  function  of  the  quadratic  integrator  with¬ 
out  Frieberger' s  pre-filter. 

Over  a  short  time  interval  (STA),  say  nAt,  an  estimate  is  made 

th 

of  the  variance.  For  the  v  epoch,  this  is  simply 
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To  take  the  place  of  the  adaptive  prefilter,  the  long  term  average 

(LTA) ,  say  m  n^t  duration.  Is  calculated  In  the  absence  of  a  seismic 

2 

event  to  provide  an  accurate  estimate  of  the  noise  variance  cr  .  The 
ratio  of  these  quantities,  namely 


Is  then  compared  to  a  constant.  For  Gaussian  random  Input  data,  this 
number  Is  the  ratio  of  two  chi -squared  random  variables  with  2 n  and 
2m  n  degrees  of  freedom.  Thus,  the  ratio  follows  the  F  distribution 
(Bendat  and  Piersol,  1971,  p.  107).  The  expected  value  of  the  ratio  Is 
unity,  and  Its  variance  Is  a  function  of  m  and  n.  For  noise  that  is  uncor 
related  with  the  seismic  signal,  the  conventional  power  detector  esti¬ 
mates  the  statistic 

S+N 

nr  * 


where  S  and  N  are  the  signal  and  noise  variances. 

If  the  STA  and  LTA  are  both  estimates  of  the  variance  of  a 
normally  distributed  quantity,  then  the  ratio  STA/LT A  is  K  -  distributed 
with  vs,vL  degrees  of  freedom  where 


vs  =  TSTA/  At 

and 

VL  =  TLTA/  At 


The  number  of  samples  in  the  STA 
The  number  of  samples  in  the  LTA 


Figure  1.4  shows  the  dependence  of  threshold  level  C  of  STA/LTA 
on  vL  for  v$  *  20  and  three  different  false  alarm  probabilities,  that  is 

Prob  [  STA/LTA  <  C]  =  PQ  . 
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STA/LTA  of  variance  estimators  for  different  false  alarm  probabilities 


Two  Interesting  features  can  Immediately  be  seen. 


1.  The  .*  is  little  to  be  gained  In  lengthening  the  LTA  beyond 
the  "knee"  of  these  curves.  For  example*  suppose  the  STA 
Is  1  second  for  20  sps  data.  Thus  vL  *  20.  With  an  LTA  of 

only  6  seconds,  vL  *  120  and  the  ratio  STA/LTA  will  be  less 
than  2.03  99  percent  of  the  time  and  less  than  2.53  99.9  per 
cent  of  the  time.  Extending  the  LTA  to  an  Infinite  period 
(perfect  estimate)  only  reduces  the  thresholds  to  1.88  and 
2.27  respectively  --  only  a  0.33db  or  0.47db  improvement. 

2.  The  change  in  the  threshold  levels  as  a  function  of  false 

alarm  probability  Is  the  other  interesting  feature.  For 
Vjj  *  20  and  »  eo  t  Table  1.1  gives  the  FAR  assuming  a 

STA  of  1.8  seconds. 


Table  1.1 


% 

CM 

1 

O 

>< 

ir> 

CM 

1 

O 

10’3 

10“4 

FAR/hr 

100 

20 

2 

.2 

Threshold 

in  db 

1.96 

2.74 

3.56 

6.5C 

Thus,  with  this  STA  period,  setting  the  threshold  at  about 
3  will  achieve  a  FAR  of  about  one  per  hour. 


The  basic  ideas  presented  here  are  all  well  known  and  form  the 
basis  of  almost  all  seismic  detectors.  The  individual  detectors  dis¬ 
cussed  in  later  chapters  differ  in  detail  and  implementation  as  each 
uses  different  methods  in  attempting  to  optimize  the  probability  of  de¬ 
tection  for  a  fixed  false  alarm  probability.  They  also  reflect  dif¬ 
ferences  in  assumed  noise  models  and  input  data  sets.  Obviously,  a 
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detector  for  low  SNR  teleselsmlc  pulses  will  be  different  from  a  high 
SNR  local  event  detector  meant  to  be  implemented  in  the  microprocessor 
of  a  field  recorder. 

1.3  TESTING  DETECTORS 

For  each  detector,  we  would  like  to  produce  the  Receiver  Operating 
Characteristic  (ROC)  curves  (Figure  1.2).  These  may  be  readily 
generated  once  the  detector  has  been  Implemented.  We  simply  generate 
a  synthetic  data  set  containing  noise  plus  signals  at  known  times 
and  with  specified  SNR.  Enough  such  data  must  be  used  to  gather  a 
meaningful  statistical  sample  of  the  detector's  performance. 

For  a  given  synthetic  test  containing  N^  events,  all  with  the 
same  SNR,  each  detector  produces  Np  false  alarms  and  NQ  detections 
of  which  Nc  are  correct.  Then 


FAR  =  Np  •  length  of  data  set/At  . 

By  repeating  the  test  with  events  covering  a  range  of  signal-to-noise 
ratios,  a  family  of  curves,  similar  to  the  single  ROC  sketched  in 
Figure  1.2,  may  be  formed. 

A  reasonable  way  to  compare  various  detection  algorithms  is  to 
test  them  against  the  theoretical  "best"  detector.  Suppose  that  the 
seismological  problems  were  completely  solved  and  the  signal  and  arrival 
time  were  perfectly  known.  Then  the  detection  algorithm  would  have  to 
choose  between 

X(t)  *  N(t)  hypotheis  Hq 

and 

X{  t)  =  S(t)  +  N(t)  hypothesis  Hj 
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where  S(t)  Is  completely  known.  Surely,  we  cannot  expect  to  do  better 
than  this. 

Helstrom  (1968)  shows  that  In  this  case,  a  monotonic  function  of 
the  likelihood  function  Is 


q(t)X(t)  dt 


(1.17) 


where  T  Is  the  detection  window  and  q(t)  Is  a  solution  of 
T 

S(t)  ■  /  q(u)  4> (t-u)  du 

J0 


(1.18) 


where  $(t-u)  is  the  covariance  of  the  stationary  gaussian  noise  N(t). 
The  expected  values  of  G  under  the  hypotheses  are 

E(G/H0)  -  0 


E(G/Hi ) 


i 

■/ 


q(t)  S(t)  dt  =  ds 


The  variance  is  the  same  under  both  hypotheses 


T  T 

Var  G  »  f  f  <l(u)  q(t)  E  [  N(u)  N(t)j  dudt 
•'q  J  0 


f  q(t)  S(t)  dt  =  d2 


(1.19) 
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by  equation  1.18  above.  The  p.d.fs  of  the  detection  statistic  are 

PQ  (G)  «  exp  [  -G2/2d2  ] 

P^  (G)  -  V2*d2  '  t  -(G-d2)2/2d2] 

and  the  false  alarm  and  detection  probabilities  are 
Q0  *  erfc  (Gg/d) 

Qq  «  1-erfc  (d-GQ/d)  (1.20) 

where  G0  is  the  decision  level  of  the  statistic  G. 

The  parameter  d  is  in  fact  the  true  signal-to-noise  ratio.  The 
easiest  way  to  see  this  is  to  consider  the  case  where  the  signal  S(t) 
is  zero  outside  the  interval  0  to  T.  Then 


OO 


S(t) 


•/ 

—  OO 


q(u)  $(t-u)  du 


Then  taking  Fourier  transforms  and  using  the  convolution  theorem 


Q(<d)  = 


and 


5(u>) 

$(u>) 


q(t) 


oO 

/ 


S(u))  elaJ^  dw 


-  OO 


*(<d) 
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Then  equation  1.19  becomes 


e1a)t  dw  I  $(t)  dt 

6*  J 


ao 


which  Is  simply  the  Integral  of  the  signal  Fourier  spectrum  divided  by  the 
noise  power  spectrum. 

Plots  of  Qd  versus  d  for  various  Q0  and  Qd  versus  Qq  for  various  d 
(the  receiver  operating  characteristics)  are  given  In  Figures  1.5  and  1.6. 
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Figure  1,5  Probability  of  detection:  completely  known  signal 
(after  Hel strom,  1978,  p.  106). 


Figure  1,6  Receiver  operatino  characteristics:  completely  known 
signal  (after  Hel strom,  1978,  p.  197). 
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II.  IMPORTANT  SEISMIC  EVENT  DETECTORS 


We  review  In  this  section  eight  seismic  detection  algorithms  which 
have  bean  developed  and  tested  by  various  organizations  over  the  past  15 
years.  Specifically,  these  are: 

1.  The  IBM  detector. 

The  Z  detector. 

3.  The  deflection  detector. 

4.  The  analytic  envelope  detector. 

5.  The  Allen  detector. 

6.  The  Stewart  detector. 

7.  The  Walsh  detector. 

8.  The  MARS  detector. 

Most  of  these  have  been  developed  for  the  detection  of  teleselsmic 
signals  recorded  by  short  period  seismometers.  Much  of  the  early  work  was 
stimulated  by  the  deployment  of  the  IASA  and  NORSAR  large  aperature  seismic 
arrays  but  more  recently  emphasis  has  been  placed  on  developing  detection 
algorithms  to  operate  on  single  traces  at  the  seismometer  location.  Such 
detectors  have  been  Implemented  on  the  SRO  seismometers  and  several  small 
portable  seismic  recording  systems.  In  these  latter  applications,  imple¬ 
mentation  Is  often  made  on  a  mocroprocessor  and  so  execution  speed  and 
algorithm  simplicity  are  at  a  premium. 

In  the  following  descriptions  of  the  detection  algorithms,  we  have 
depended  heavily  on  the  reference  given  at  the  beginning  of  each  sub¬ 
section.  Unless  otherwise  stated,  text  enclosed  in  quotes  is  taken  from 
these  reports.  In  the  references  section  at  the  end  of  Section  IV,  all 
quoted  literature  is  listed  along  with  the  abstracts  of  the  principal 
works.  In  Appendix  A,  a  more  complete  bibliography  is  to  be  found. 

2.1  IBM  DETECTOR,  VANDERKUIK,  et  al_.  (1965) 

This  report  describes  the  seismic  detv  .tion  algorithm  developed  in 
1965  by  IBM  for  use  on  the  LASA  array  and  subsequently  implemented  at 
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NORSAR.  The  LASA  array  consisted  of  525  seismometers  grouped  Into  21  sub¬ 
arrays  of  25  seismometers  each.  The  computer  system  used  for  detection 
has  as  Its  Input  some  300  beam  outputs,  the  subarray  vertical  beams,  and 
the  outputs  from  21  single  seismometers  (one  from  each  subarray)  attenuated 
to  record  large  signals. 

"The  processing  applied  to  a  single  LASA  beam  output  can  be 
described  as  a  filtering  process  designed  to  de-emphaslze  those  frequency 
components  where  the  beam  output  slgnal-to-nolse  ratio  Is  low  compared  to 
Its  maximum.  The  filtering  process  Is  followed  by  an  Integration  opera¬ 
tion  In  which  either  the  square  or  the  magnitude  of  the  filtered  beam 
channel  Is  Integrated  over  a  time  interval  of  fixed  duration.  The  re¬ 
sulting  single  quantity  must  then  be  compared  with  a  threshold  value  for 
detection  purposes." 


f i 1  ter 


i 


l 


recti fy 


Integrate 


Figure  2.1.  The  IBM  detector. 


>C 

<C 


This  detector,  schematically  represented  by  a  flow  diagram  In  Figure 
2.1,  Is  quite  close  to  the  one  prescribed  by  Frieberger  (1963)  as  described 
in  Section  1.2,  except  that: 

1.  The  filter  used  was  not  the  optimum  filter, 


I 

V  fmTOCFsW 

but  simply  ^"T7f^(w)  which  is  the  large  SNR  limit  of 
the  optimum;  (In  fact,  according  to  Blandford  (1980  Private 
Communication),  IBM  did  not  use  an  adaptive  filter  or  even 
a  1/Fn(w)  filter  but  merely  a  bandpass.) 

2.  The  absolute  value  of  the  filter's  output  was  taken  for  compu¬ 
tational  efficiency  (speed)  rather  than  the  square; 
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3.  The  integrator  was  a  "leaky"  integrator  Implemented  by  a 
digital  recursive  filter  rather  than  a  proper  Integrator. 

Vanderkulk  (1964)  examined  the  effects  of  these  three  approximations. 
The  filter  with  an  amplitude  response  of  l//FN(al)  was  applied  only  over 
a  band  <Dj  <  w  _<  oc^  and  zero  outside.  The  phase  was  arbitrary.  To  gauge 
the  effects  of  such  a  filter,  the  report  assumed  that  Fs(w)/FN(w) 
paaked  at  frequencies  In  the  center  of  the  band  and  fell  off  ex¬ 

ponentially  on  either  side  so  the  ratio  was  L  db  below  the  peak  at 
and  uX). 

"Figure  2,2  depicts  the  graph  of  the  loss  versus  L.  As  could  be 
expected,  the  graph  shows  a  minimum:  when  L  is  small,  the  performance 
of  the  noise-prewhltenlng  filter  suffers  because  too  little  of  the 
signal  is  passed;  when  L,  and  hence  the  bandwidth,  is  large,  the  per¬ 
formance  deteriorates  because  the  filter  passes  too  much  noise. 

Note  that  the  loss  in  performance  Is  less  than  1  db  when  L  Is 
between  2.2  db  and  12  db.  Thus,  to  assure  a  loss  of  less  than  1  db,  the 
nolse-prewhitenlng  filter  bandwidth  must  be  large  enough  to  include  all 
frequencies  where  the  input  signal-to-noise  ratio  is  within  2.2  db  of 
its  maximum,  but  must  reject  all  frequencies  where  the  signal-to-noise 
ratio  is  more  than  12  db  below  the  maximum.  Thus,  in  practice,  a 
comfortable  frequency  margin  is  allowed  in  which  to  achieve  the  filter 
cutoff.  The  above  result  permits  another  interpretation  as  well.  It 
shows  that  those  frequency  components  of  the  input  channel  where  the 
signal-to-noise  ratio  is  more  than  2.2  db  below  the  peak  value  are 
ineffective  so  far  as  signal  detection  is  concerned.  Therefore,  the 
signal  processing  which  produced  the  input  channel  is  allowed  to  be 
arbitrarily  degraded  for  those  frequencies  where  the  signal-to-noise 
ratio  is  more  than  2.2  db  below  the  maximum,  so  long  as  the  signal-to- 
noise  ratio  in  the  region  where  it  is  within  2.2  db  of  the  maximum  is 
not  materially  altered." 

The  effect  of  rectifying  rather  than  squaring  the  output  of  the 
prefilter  was  also  investigated  in  this  report. 
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Figure  2.2.  Loss  In  performance  of  the  noise-prewhitening  filter 
for  the  IBM  detector  (after  Vanderkulk,  et  al.,  1965). 
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"When  rectifying  and  integrating,  the  output  produced  is  given  by 


dt  . 


(2.1) 


When  squaring  and  integrating,  the  output  is  given  by 


T 

P  •  }  f  x(t)2  dt  .  (2-2) 

J0 

The  channel  trace,  x(t),  Is  assumed  to  be  the  sum  of  Gaussian  stationary 
signal  (when  present)  and  Gaussian  stationary  noise.  The  loss  incurred 
by  applying  rectified  integration  instead  of  squaring  and  integrating 
is  found  by  dividing  the  signal-to-noise  ratio  of  P  by  that  of  R."... 
"The  resulting  quotient,  Q,  is  given  by 


Q  -  i±4zhm  A , 


(2.3) 


.here 


T  T 


A2  = 


tt  -  2)  f  f  p(t2  -  t,)  dt1  dt2 

Aro _ 

T  T 

f  f  p(t2  ■  t-j )  dt.j  dt2 

J  n  J 


(2.4) 


in  which  formula  p(t2  -  t-)  designates  the  correlation  coefficient  be- 
tween  x ( t^ )  and  x(t2)  when  noise  alone  is  present,  and  p(t2  -  tj) 
designates  the  correlation  coefficient  between  |x(t^)|  and  |x(t2)| 
when  noise  alone  is  present.  Furthermore,  S  is  the  total  signal 
power,  i.e., 


S  =  f  S(f)  df  . 
J0 


'2.5) 


20 

S YSTSMM.  3CICHGK  AS  '  SOPTWAt te 


Likewise,  N  Is  the  total  noise  power.  Thus,  S/N  may  be  termed  the 
Input  slgnal-to-nolse  power  ratio. 

It  Is  readily  shown  that 

p  *  arcsine  p  +^1  -  p2  -  lj/(j  -  l)  ,  (2.6) 

from  tah-ich  It  follows  that 

p2/(t  -  2)  <  a  ip2  .  (2.7) 

Consequently,  1  <  A  <  A  -  2  =1.07.  Since  the  loss  in  performance  Is 
10  1o9-|qQ  db,  this  loss  is  given  by  the  following  expression: 

Loss  =  10  log1Q  p-  xl.035j  db  ,  (2.8) 

where  the  error  due  to  equating  A  to  1.035  is  less  than  0.2  db."  For 
this  model  the  loss  in  performance  is  not  significantly  affected  by  the 
shape  of  the  noise  spectrum  or  by  the  choice  of  the  intetration  time  T. 
A  graph  of  the  loss  versus  the  input  signal-to-noise  ratio  S/N  is 
shown  in  Figure  2.3. 

Finally,  the  effects  of  the  leaky  integrator  were  examined. 

Most  of  the  detectors  that  are  implemented  using  sampled  data  use 
recursive  filters  to  simulate  integrators.  Typically,  the  output 
at  time  t^  is  given  by 

Y.  *  (1  -  OY.^  +  C  X.  (2.9) 

where  C  is  a  constant  <1  and  X.  is  the  input.  In  matrix  notation 
this  car)  be  approximated  by 

Y  =  a  M  X  (2.10) 
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1 


where 


1 

0 

0 

0 

0  • » • 

(1  -  C) 

0 

0 

0 

o  •  •  * 

H  * 

(1  -  c)2 

(1 

- 

C) 

1 

0 

o  . .. 

(2.11) 

(1  -  C)3 

• 

(1 

• 

C)2 

(1  -  C) 

• 

1 

• 

o  •  •  • 

• 

j 

• 

• 

• 

• 

• 

• 

• 

• 

•  •  •  • 

• 

The  exponential  integrating  filter  output  can  be  represented  by 
t 

Y  =  J  x(u)2  exp  (-  t  j  u)  du  , 

•  00 

while  the  straightforward  power  integration  is  given  by 


(2.12) 


t+T 

1=1  x(u)2  du  . 

't 


■/ 


(2.13) 


In  order  to  assess  the  effect  of  finite  signal  duration,  the  signal 
is  taken  to  be  a  portion  of  duration  of  a  Gaussian  stationary  process. 
It  follows  that  optimum  processing  (maximized  signal-to-noise  ratio)  is  pro¬ 
duced  by  using  the  output  Z  with  T  =  Ts  and  with  the  interval  from  t 
to  t  +  T  exactly  covering  the  signal  time  interval. 

With  this  maximum  possible  signal-to-noise  ratio  as  a  standard, 
the  loss  in  performance  resulting  from  the  use  of  the  exponentially  de¬ 
caying  integrating  filter  (with  T  being  arbitrary)  is  given  by  the 
following  expression: 


Loss  =  -10  log1Q  (l  -  exp(-T$/T)) 


(2.14) 
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This  formula  was  obtained  by  computing  the  quotient  of  the  maximum 
possible  signal -to-noise  ratio  and  the  signal-to-noise  ratio  of  the 
rectified  integration  output  Y,  with  t  being  the  endpoint  of  the 
siynal  time  interval.  The  loss  is  than  equal  to  10  times  the  logarithm 
of  this  quotient.  The  preceding  formula  for  the  loss  is  an  approximation 
obtained  by  assuming  that  both  T$  and  T  are  large  compared  with  the 
reciprocal  of  the  noise  bandwidth  (or,  equivalently,  compared  with  thp 
decay  time  of  the  noise  autocorrelation  function  p ( t) ) .  Figure  2.4 
is  a  graph  of  the  loss  versus  T/Ts.  The  minimum  loss  occurs  when 
T  *  0.80  Ts  and  is  0.45  db.  The  performance  loss  is  less  than  0.55  db 
when  T  lies  anywhere  between  0.55  Tg  and  1.10  T  .  Thus,  by 
tolerating  this  loss,  it  is  sufficient  to  implement  this  exponentially 
decaying,  integrating  filter  with  a  sequence  of  filter  decay  times  T 
which  progresses  by  factors  of  two. 

The  separate  performance  loss  considerations  given  above  should 
more  properly  have  been  combined  to  provide  the  simultaneous  effect 
of  these  three  factors  on  the  signal  processing  performance. 

The  only  difference  between  the  "IBM"  detector  and  the  conven¬ 
tional  power  detector  (see  Section  1.2)  is  that  the  former  takes  the  posi- 
tive  square  root  of  (i.e.,  the  absolute  value,  a  computationally 
efficient  quantity)  rather  than  the  square.  As  we  saw  earlier,  under  a 
simple  model  of  the  noise  signal,  the  loss  in  performance  caused  by 
taking  | |  instead  of  |x^|  is  only  about  1  db  (0.05  mb  difference  in 
detection  capability). 

The  IBM  approximation  to  the  conventional  power  detector  is 
basically  the  algorithm  that  is  used  at  both  N0RSAR  and  SDAC. 

The  data  sampled  at  20  sps  is  first  decimated  (without  filtering) 
to  10  sps.  The  data  is  next  passed  through  a  recursive  band  pass  filter 
and  further  decimated  to  5  sps.  This  is  the  input  data  stream  to  the 
rectifier  (absolute  value) 

=  |x .  |  where  Ai  =0.2  sec.  (2.15) 
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Loss  (db) 


these  rectified  quantities  are  then  averaged  Into  groups  of  three 


1 

where  now  Ai  corresponds  to  (2.16) 
0.6  sec,  l'.e. ,  a  decimation 
by  3. 

The  three-point  averages  are  then  further  averaged  over  three 
0.6  second  Intervals  (1.8  second)  (i.e.,  no  decimation)  to  form  the 
STA,  which  Is  updated  every  0.6  second. 


1 


j=1-2 


Aj  corresponds  to  0.6  sec 
Ai  corresponds  to  0.6  sec 


(2.17) 


i 


The  long-term  average,  LTA,  is  formed  by  a  first-order  recursive 
filter  acting  on  the  0.6  second  samples  of  the  STA.  However,  there  is  a 
further  decimation  by  three  so  that  the  LTA  filter  output  is  derived 
from  statistically  independent  samples  of  the  STA 


LTA 


.  =  (1-C)  LTA._1  +  C  STA 


3j 


Aj  corresponds  to  0.6  sec 
Ai  corresponds  to  1.8  sec 


Typically,  C  is  set  for  a  30  second  time  constant: 


c  *  ff  *  ir =  °-06- 


(2.18) 


(2.19) 


The  decision  algorithm  examines  the  0.6  second  samples  of  STA 
and  compares  this  to  the  1.8  second  samples  of  the  LTA.  If 

STA  >  K  *  LTA  (2.20) 

for  Q  out  of  Q'  successive  tests  an  event  is  declared.  Typically 
Q/Q'  is  set  at  2/3  or  3/3  and  K  is  approximately  3.  Note  that  since 
LTA  is  only  updated  every  1.8  second,  it  remains  constant  over  this 


test. 

When  an  event  is  declared,  the  time  constant  of  the  LTA  re¬ 
cursive  filter  is  reduced  typically  to  half  its  normal  value  until  the 
event  is  declared  over.  This  is  done  to  detect  the  later  phases. 
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2.2 


"Z"  DETECTOR ,  SWINDELL  AND  SNELL  (1977) 


The  "Z"  detector  was  the  result  of  research  to  develop  a  constant 
FAR  (CFAR)  algorithm  for  the  Station  Processor.  This  "detector  is  a 
modification  of  a  conventional  power  detector  which  detects  signals  as 
short-term-averaged  (STA)  signal  power  relative  to  long-torm-averaged 
(LTA)  noise  power  preceding  the  signal.  Statistical  analyses  of  STA 
noise  fluctuations  indicated  that  small  deviations  in  the  STA  standard 
deviation  from  LTA  causes  serious  problems  In  controlling  false  alarms. 
This  sometimes  leads  to  the  unstable  operation  of  convenl tonal  power  de¬ 
tectors.  In  some  cases,  such  instability  causes  the  detector  to  turn  on 
or  to  shut  off  for  long  periods  of  time.  The  "Z"  detector  was  designed 
to  solve  this  problem  by  continuously  adjusting  the  threshold  of  STA- 
LTA  to  a  fixed  number  of  standard  deviations  of  STA-LTA.  Estimates  of 
the  standard  deviation  are  updated  on  a  point-by-point  basis.  The  "Z" 
detector  is  also  designed  to  control  false  alarms  from  the  coda  of  large 
'ignals."  Basically,  the  detector  statistic  for  noise  input  was  modeled 
s  a  log-normal  process.  The  mean  and  variance  of  the  statistic  were  used 
>  transform  the  statistic  to  a  zero-mean,  unit  variance  quantity  called 
the  "Z"  statistic.  Thus,  for  a  single  random  variable  x,  the  Z 
statistic  is  defined  as 

(2-a) 

The  idea  for  the  logarithmic  transformation  came  originally  from  LaCoss's 
(19:2)  observation  that  the  STA  values  from  the  IBM  detector  looked  to 
be  normal  distributed.  Of  course,  if  the  input  signal  were  Gaussian, 

and  the  STA  values  were  a  true  estimate  of  the  power,  they  would  be 

2 

X  distributed.  This  distribution,  for  a  small  number  of  degrees 
of  freedom,  is  not  unlike  a  log  normal  distribution  in  appearance. 

However,  as  the  IBM  STA  algorithm  does  not,  in  fact,  estimate  power,  cut 
rather  averages  the  absolute  value,  LaCoss's  empirical  approach  to  the 
correct  distribution  is  reasonable. 

In  any  event,  the  "Z"  detector  proceeds  as  follows: 

1.  The  input  data  are  band  passed  with  a  recursive  filter. 

This  filter  is  centered  at  about  2  Hz  and  is  between 
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3  and  4  Hz  wide.  The  principal  purpose  of  the  filter 
Is  to  remove  the  mean  from  the  data.  It  is  not  chosen 
"a  la  Fr1eberger"to  optimize  detectibility. 

2.  The  beginning  operation  for  the  detection  statistic  Is 
a  squaring  rather  than  rectifying  operation 

Y1  »  .  (2.22) 


3.  As  In  the  NORSAR/SPAC  Implementation  of  the  conventional 
power  detector,  the  Y.  are  averaged  over  some  gate  N 
by 


j 


1*j-N 

The  time  interval  associated  with  Y  is  N  times  longer 
than  that  associated  with  Y.  M  sets  of  these  points 
are  then  averaged  to  form  the  STA  which  yields 


(2.24) 


STA  is  updated  just  as  frequently  as  Y. 

4.  "The  logarithm  of  the  short-term  average,  STA,  is  passed  to 
the  long-term  background  level  estimator  to  be  used  for  up¬ 
dating  the  long-term  estimate,  u,  of  log  STA  and  of  the  vari- 
2 

ance  a  ,  of  log  STA.  Z  is  computed  using  p  and  a  before 
they  are  updated  with  the  present  value  of  log  STA.  Then  a 
sequential  threshold  test  is  performed.  The  first  test  com¬ 
pares  the  absolute  magnitude  of  z  to  the  threshold  value 
of  z,  ZyH.  This  is  to  prevent  the  LTA  from  updating  either 
on  true  detections  (positive  z)  or  highly  negative  values 
of  z  arising  from  the  logarithm  of  very  low  powers  which 
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occur  momentarily.  If  the  threshold  Is  exceeded,  the  sign 
of  z  Is  checked  and,  If  positive,  a  detection  Is  declared." 

5.  The  fifth  step  Is  a  detection  analyzer  which  "suppresses  multiple 
declarations  of  detections  from  the  same  signal  caused  by  coda 
levels  continuously  exceeding  the  detection  threshold,  posts 
new  detection  declarations  If  coda  levels  show  unusual  in¬ 
creases  In  level,  and  controls  the  IT A  estimator  to  prevent 
LTA  estimates  from  being  contaminated  by  coda  energy." 

The  concept  of  a  post  detection  processor,  found  In  other  detectors 
as  well.  Is  Important.  It  allows  one  to  build  a  "front  end"  system  that 
is  fast,  which  has  a  large  FAR,  but  which  culls  the  data  stream  to  a 
more  manageable  volume.  In  the  "Z"  detector,  "the  internal  logic  of  the 
analyzer  Is  moderately  complicated  and  Is  best  understood  by  Its  action 
on  a  typical  signal  (see  Figure  2.5).  It  may  be  broken  down  further  into 
a  coda  suppressor  and  LTA  controller. 

When  the  input  z  crosses  the  threshold  zTH,  the  time  is  noted 
and  a  timer  is  initiated.  After  the  specified  time,  the  beam  select 
gate,  has  elapsed,  the  detection  time  and  the  peak  value  of  z  occurring 
in  this  gate  is  stored  for  the  beam  selector's  use.  The  beam  selector 
gate  is  operator  adjustable  and  typically  lasts  ten  seconds.  Also  at  the 
detection  time,  a  second  timer  is  Initiated  which  defines  a  secondary 
detection  gate.  (This  gate  is  also  operator  adjustable;  during  this  study, 
15  seconds  was  used.) 

Assume  for  the  moment  that  the  signal  is  small  and  the  coda 
level  drops  below  the  threshold  a  few  seconds  after  the  initial  detection. 
The  situation  now  is:  no  detections  are  occurring  at  the  z  comparator, 
there  is  remnant  signal  coda  energy  present,  but  it  is  decaying  toward 
the  original  noise  level.  When  z  becomes  less  than  z^,  a  third  timer, 
also  operator  adjustable,  is  initiated  which  defines  the  detector  reset 
time.  After  this  time  has  elapsed,  the  detector  is  reset  and  any  new 
threshold  crossing  will  be  declared  a  new  detection.  For  this  study,  a 
reset  time  of  20  seconds  was  used.  The  purpose  of  the  reset  time  is  to 
prevent  new  detections  from  being  declared  because  of  "jitter"  in  the 
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.  Coda  supressor  for  the  "Z"  detector  (aftosr  Swindell  and 
Snell,  1977). 
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z  values  as  they  decay  through  the  threshold.  The  reset  time  must  be 
greater  than  the  secondary  detection  gate  which  In  turn  Is  greater  than 
the  beam  selector  gate. 

The  reaction  to  a  large  signal  where  z  remains  above  z^  for 
at  least  two  secondary  detection  gate  Intervals  Is  slightly  more 
elaborate.  The  goal  is  to  monitor  the  coda  level  for  the  arrival  of 
secondary  phases  or  signals  from  another  event.  The  operating  assumption 
Is  that  the  coda  level  will  decay  fairly  smoothly  from  its  Initial  peak 
except  when  new  signals  arrive.  As  each  secondary  detection  gate  lapses 
(except  the  first  one),  the  peak  z  in  that  gate  is  compared  with  the 
peak  value  from  the  preceding  gate.  If  it  exceeds  the  old  peak  by  some 
specified  amount  (e.g.,  6  dB,  operator  adjustable)  than  a  new  detection 
is  declared. 

The  LTA  controller  maintains  the  highest  level  of  control  of  the 
LTA  estimator  and  its  job  is  to  exclude  as  much  signal -related  energy  as 
possible  from  the  estimates  of  the  mean  and  standard  deviation  of  the 
log  STA  for  the  background  noise.  Its  action  can  also  be  described  more 
easily  by  example  (refer  to  Figures  2.6  and  2.7).  After  computing  z, 
the  LTA  estimator  normally  updates  u  and  a2  with  every  new  datum 
except  when  inhibited.  The  threshold  comparator  issues  an  inhibit 
command  whenever  |z|  >  z^.  Finally,  to  insure  that  a  very  large  signal 
does  not  keep  y  frozen  too  long  the  freeze  time  is  limited  to  some 
maximum,  say  10  minutes,  which  is  also  operator  adjustable.  When  y 
becomes  unfrozen  by  any  means,  it  assumes  the  present  value  of  LTA 
and  normal  updating  resumes.  Internal  to  the  LTA  estimator,  there  are 
actually  two  separate  estimates  of  mean  and  standard  deviation  of  log  STA. 
Using  the  mean  as  an  example  (the  o2  estimate  follows  in  parallel), 
the  two  quantities  are  y  and  "LTA".  Z  is  always  computed  using  y. 

Under  conditions  of  noise  with  no  inhibiting  commands  from  the 
LTA  controller,  y  and  LTA  are  identical.  When  a  detection  occurs,  the 
dichotomy  of  LTA  and  y  becomes  apparent.  Whenever  |z|  <  ZyH,  LTA 
updates  but  y  remains  frozen  at  least  until  the  L^A  reset  time  interval 
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Figure  2.6.  "Z"  detector  flow  diagram  (after  Swindell  and 

Snell,  1977). 
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Figure  2.7.  "Z"  detector  LTA  controller  (after  Swindell 

and  Snell ,  1977). 
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has  elapsed.  At  the  end  of  the  reset  time  a  test  of  the  difference  be¬ 
tween  LTA  and  p  is  made  and  y  remains  frozen  until  that  difference 
becomes  small  enough.  This  maximum  difference  is  called  the  secondary 
reset  function  limit  and  is  operator  adjustable." 

Evaluation  of  the  "Z"  detector  was  reported  in  Swindell  and  '„nell 
(1977)  using  data  from  the  KSRS  array  during  different  periods  a  the 
year.  They  concluded  that  "the  Z  statistic  algorithm  produces  a 
constant  false  alarm  rate  rather  than  a  constant  alarm  Kute  and  Is  es¬ 
sentially  Independent  of  the  noise  field  behavior,  "r,ie  estimated  mean 
and  variance  of  the  noise  which  is  used  to  com/1-,  c  the  basic  detector 
output  to  Z  also  provides  information  u*-.,ul  for  estimating  station 
performance  :r  *  — «na  makes  unnecessary  a  separate  noise 

level  calculation  for  the  power  detector.  Alarm  rates  may  be  set 
independently  for  each  beam  (signal)  to  reflect  differences  in  seismicity 
or  to  give  some  beams  (signals)  higher  sensitivity  without  affecting 
overall  alarm  rate." 

They  also  tested  two  different  pre-detection  band  pass  filters 
and  concluded  that  they  caused  no  appreciable  differences  on  the  out¬ 
come,  nor  did  changes  in  the  integration  time  of  the  STA  process  from  1.6 
second  to  3.2  seconds  have  appreciable  effect. 

Further  tests  of  this  automatic  detector  were  conducted  during 
two  months  of  field  operation  at  DET  459  by  Secoy  (1978).  Sax  (1980) 
reports  that  "at  KSRS,  70  percent  of  the  analyst  picks  were  automatically 
detected  with  0.875  false  alarms  per  hour  and  with  a  stability  of  about 
35  percent.  At  DET  459,  70  per  cent  of  the  analyst  picks  were  automatically 
detected  with  0.375  false  alarms  per  hour.  These  false  alarm  figures  are 
on  a  per-beam  basis  so  that  they  are  comparable  with  those  expected  for 
a  single-sensor  channel.  The  timing  accuracy  of  the  detector  indicated 
a  late  bias  of  1.7  seconds,  with  a  standard  deviation  of  2.2  seconds. 

The  potential  exists  for  making  major  improvements  in  the  "Z" 
detector.  A  need  exists  to  reduce  false  alarms  caused  by  local  seismic 
signals,  and  to  provide  a  restart  mechanism  for  automatically  "warming- 
up"  the  detector  after  long  power  outages.  Also,  there  is  a  need  to 
provide  post-detection  processing  in  order  to  improve  the  timing  accuracy 
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of  signals,  and  possibly  to  cope  with  multiple  arrivals  of  complex  and 
multlpathed  signal  transmissions.  A  recent  application  of  the  "Z"  de¬ 
tector  by  Sax,  et  al_.  (1979b)  Indicates  that  the  detector  at  KSRS  pro¬ 
vided  precise  magnitude  estimates  of  complex  events  from  a  selected 
USSR  border  region.  Timing  associated  with  event  arrivals,  In  this 
case,  was  complicated  by  the  repeated  arrival  of  other  events,  at 
equal  perceived  magnitudes,  for  up  to  a  minute  following  the  first  ar¬ 
rival  . 

Additional  work  Is  needed  to  develop  pre-filters  to  optimize  the 
detection  of  earthquakes  and  explosions.  Since  the  signal-to-nolse 
ratio  (SNR)  for  explosions  and  local  events  may  be  significantly  different 
from  that  for  earthquakes,  it  may  be  desirable  to  provide  more  than  one 
optimized  pre-filter  for  more  complete  extraction  and  identification  of 
various  signal  types. 
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2.3  DEFLECTION  DETi-CTOR,  SHENSA  (1977) 


The  deflection  detector  works  In  the  Fourier  domain  by  taking 
data  vectors  x  and  calculating  the  digital  Fourier  transform  of 
that  data  set.  Using  an  FFT  algorithm,  the  power  P^(k)  (for  each 
frequency  k,  P. (k)  Is  the  real  part  of  the  FFT  squared  plus  the 
imaginary  part  of  the  FFT  squared)  at  the  k  frequency  is  found  for 

the  ith  epoch  or  time  wondow.  This  moving  power  spectrum  estimator  is 
used  in  three  closely  related  detectors  (see  Figure  2.8): 

1.  The  average  (or  stacked)  power  detector. 

2.  The  maximum  deflector  detector. 

3.  The  average  (or  stacked)  deflection  detector. 

In  each  of  these  cases,  Shensa  (1977)  used  both  P^(k)  and  P."(k)  = 
log  Pj(k)  as  the  basic  input  quantity,  but  in  all  cases  he  found  that 
taking  the  logarithm  decreased  the  probability  of  detection  for  a  given 
SNR. 


2.3.1  The  Average  Power  Detector 

This  detector  is  simply  the  "Z"  detector  operating  in  the  frequency 
domain.  The  average  of  P^(k)  over  some  signal  band  is  found  and  a  Z- 
statistic  calculated  by  normalizing  to  the  mean  and  standard  deviation 
of  the  average  power 

N2 

*E  pi<k>  -» 

k*N, 

Y.  =  - L -  N  =  N2  -  Nx  (2.25) 

where  p  and  o  are  the  mean  and  standard  deviation  of  the  average 
power. 

This  is  really  an  estimate  of 
Var.  -  u(var) 

o(»r)  •  <2-26> 
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Figure  2.8.  Block  diagram  of  three  deflection  detectors  showing  how  the 
scalar  statistic  X,  Y  or  Z  is  derived  from  the  power  spectrum 
Pj(k)  (After  Shensa,  1977). 
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.arxTT.-^-jB &/zxm:xr.  2WBB  aCB 


s.'  ■  r?T,\  *i 


or  simply  a 


STA(var)  -  LTA(var) 
a(var) 


(2.27) 


which  Is  just  what  Swindell  and  Snell  (i.277)  calculated  in  the  time 
domain  for  the  "Z"  detector,  after  taking  the  log  of  the  short  term 
variance  estimator.  An  event  is  declared  if  this  ratio  exceeds  some 
fixed  threshold. 


2.3.2  The  Maximum  Deflector  Detector 

The  power  (k)  is  converted  into  a  "Z- statistic"  on  a  fre¬ 
quency-  by-  frequency  basis  by  the  usual  transformation 


Zi(k) 


Pj(k)  -  ti(k) 


(2.28) 


where  y(k)  and  a(k)  are  the  mean  and  standard  deviation  of  the  power 
at  the  kth  frequency  calculated  in  the  absence  of  a  seismic  signal.  At 
each  epoch  (value  of  i)  the  maximum  Z..(k)  is  compared  against  a  fixed 
threshold  value  and  an  event  declared  if  it  exceeds  this  value. 


2.3.3  The  Average  Deflection  Detector 


In  this  detector,  instead  of  choosing  the  max  [Z^ (k)]  the 
average  across  some  band  is  used 


Xi  =  FT  Zj(k)  where  N  *  (2.29) 

i=Ni 

With  this  implementation,  the  standardized  variable  Z^  is  computed 
on  a  frequency  band-by- frequency  band  basis.  There  is  a  certain  simpli¬ 
fication  by  noting  that  for  the  unsmoothed  spectral  estimators  used  here, 
the  standard  theory  of  random  data  analysis  yields  the  result  that  in 
the  absence  of  signal  the  sampling  distribution  of  P..(k)  is  given  by 

2  Mk) 

P^k)  =  X|  -  (2-30) 
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where  is  the  chi  squared  random  variable  with  two  degrees  of  freedom, 
and  (k)  is  the  true  power  spectral  density.  Since 


we  see 


y  ~  P  ^  ( k )  > 
a  ~  P^(k)  . 


This  is  the  well  known  result  that  in  the  absence  of  smoothing,  the 
standard  derivation  of  a  power  estimate  equals  the  estimate  itself. 
Thus,  equating  the  LTA  with  y  for  each  frequency,  we  get 


(2'31) 

Although  P.(k)  is  itself  distributed,  the  statistic  on  Z.(k) 

is  not  so  simple  because,  of  course,  y  and  o  themselves  are  not 
known  exactly  but  must  be  estimated  from  the  data. 

This  detector  statistic  is  proportional  to  a  weighted  power  spec¬ 
tral  density  average  between  N-j  and  It,.  In  contrast  with  the  average 
power  detector,  the  average  deflection  detector  weights  each  power 
spectral  density  component  inversely  with  the  estimated  standard  devia¬ 
tion  of  the  noise  power  spectral  density.  The  detection  statistic  used 
produces  the  same  incoherent  noise  gain  as  the  average  power  detector, 
but  it  is  optimum  only  for  indpendent  spectral  noise  fluctuations.  The 
average  deflection  detector  differs  significantly  from  the  average 
power  detector  in  that  it  weights  the  signal  spectrum  more  heavily  at 
those  frequencies  where  it  peaks  relative  to  the  noise  spectrum.  Note 
the  similarity  between  this  and  the  ^  l/FN(uj  prefilter  of  the  IBM 
detector  (Section  2.1),  and  the  prewhitening  operation  of  the  Walsh 
detector  (Section  2.7). 
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Several  weak-signal  scenarios  are  given  to  illustrate  factors  which 
influence  the  comparative  performance  of  the  above  three  detector  concepts. 

The  first  hypothesis  we  shall  consider  is  that  weak  signals  contain 
energy  which  significantly  exceeds  the  noise  in  at  least  one  of  the  fre¬ 
quency  bands  monitored  by  the  detector.  In  this  case,  the  maximum  deflec¬ 
tion  detector  will  efficiently  detect  such  a  coherent  sinusoidal  signal. 

The  average  deflection  and  average  power  detector  would  more  than  likely 
miss  this  type  of  signal  since  averaging  over  other  frequency  cells 
would  reduce  this  signal  peak  by  as  much  as  1/N  (N  frequency  cells) 
while  the  incoherent  gain  in  reducing  the  noise  standard  deviation  would 
be  no  more  than  /fT.  Thus  a  detection  loss  as  large  as  1/.4T  from  the 
coherent  gain  of  the  maximum  deflection  would  result  from  averaging  all 
frequencies  over  the  full  band  of  the  detector. 

A  second  hypothesis  is  that  the  weak  signal  spectrum  does  not  exceed 
the  noise  by  a  significant  amount  in  any  one  of  the  frequency  bands,  but  it 
is  near  or  slightly  above  the  noise  spectrum  in  the  entire  band  covered 
by  the  detector.  In  this  case,  the  maximum  deflection  detector  would  miss 
the  signal;  in  addition,  averaging  over  the  band  would  produce  a  /N 
incoherent  gain  by  the  average  power  detector.  Averaging  would  produce 
a  somewhat  smaller  gain  in  the  average  deflection  detector  in  that  a 
statistically  variable  weighted  sum  would  degrade  the  detector  perform¬ 
ance.  This  expected  incoherent  gain  of  VTT  would  be  produced  by  re¬ 
ducing  the  standard  deviation  of  noise  fluctuations.  The  latter  ;ould 
be  accomplished  by  averaging  the  power  spectral  densities  of  the  N 
frequency  components  covering  the  signal  band. 

The  third  hypothesis  is  similar  to  the  second  in  that  the  power 
spectral  density  of  the  signal  is  near  that  of  the  noise.  Consequently, 
no  detection  is  possible  by  the  maximum  deflection  detector  because  no 
large  coherent  gain  can  be  observed  at  any  frequency.  However,  in  this 
case,  the  noise  spectrum  is  unstable,  and  it  varies  greatly  with  statis¬ 
tical  independence  from  one  frequency  cell  to  another,  and  also  from  one 
time  slice  to  another.  Similarly,  the  signal  power  fluctuations  might 
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vJifc, 


also  be  unstable.  In  either  case,  the  weighted  average  of  the  stacked 
deflection  detector  would  be  required  In  order  to  achieve  the  most  effi¬ 
cient  detection  of  weak  signals. 

The  conditions  for  optimum  performance  of  the  three  types  of 
spectral  detectors  are  summarized  as  follows: 

t  The  maximum  deflection  detector  Is  optimized  for  a  wcJ<  signal 
If  the  signal  significantly  exceeds  the  noise  at  only  one 
or  several  frequencies  within  the  signal  band. 

•  The  average  power  detector  is  an  optimum  detector  jf  the 
weak-signal  and  noise  spectra  are  stable  and  if  the  signal 
spectrum  is  near  or  slightly  exceeds  the  noise  spectrum 
uniformly  over  a  broad  range  of  frequencies. 

•  The  average  deflection  detector  is  an  optimum  detector  if 
the  weak-signal  and  noise  spectra  are  highly  unstable  and 
if  the  signal  spectra  is  near  or  slightly  exceeds  the  noise 
spectrum  over  a  broad  range  of  frequencies. 

Shi  a  evaluated  the  performance  of  these  three  detectors  by  adding 
■'.wo  earthquake  samples  and  two  explosion  samples  to  3,000  seismic  noise 
samples.  He  examined  the  detection  performance  at  various  input  SNRs 
(SNRs  between  -6  and  +9  dB).  His  goal  was  to  obtain  the  complete  operating 
characteristics  (P^  versus  Pf)  for  all  three  types  of  detectors.  The 
most  comprehensible  results  were  obtained  by  comparing  the  detectors  for 
32-point  FFTs  and  P^-versus-SNR  for  a  false  alarm  rate  constrained  at 
4.5  false  *  arms  t  hour. 

In  the  case  of  all  four  signals,  the  average  deflection  detector 
operating  characteristics  are  much  worse  than  are  either  the  average 
power  or  the  maxirn  'eflection  detector.  This  indicates  that  the  seismic 
signal  and  noise  ,  .otra  are  nearly  stable  and  stationary.  These  results 
also  indicate  that  the  average  of  spectral  "Z"  statistics  is  neither 
optimum  nor  feasible  for  use  as  a  detector  statistic. 


41 


SYSTEMS.  SCIENCE  AND  SOFTWARE 


In  the  case  of  the  two  earthquake  signals,  at  all  input  SNRs,  the 
average  power  detector  performed  best.  At  4.5  false  alarms  per  hour  and 
for  a  detection  /robability  of  0.5,  the  average  power  detector  gained 
0.2  mb  of  detection  capability  over  the  maximum  deflection  detector 
in  one  case,  and  by  0.05  mb  in  the  other  case.  Thus,  the  broadband  "Z" 
detector  appears  to  be  a  superior  detector  of  earthquakes  but  not 
necessarily  of  explosions.  These  operating  results  are  relevant  to 
detecting  weak  signals  (that  is,  signals  which  are  only  about  6  dB 
above  the  ambient  noise). 

In  the  case  of  the  two  explosion  signals  and  at  all  input  signal  - 
to-noise  ratios,  the  maximum  deflection  detector  performed  significantly 
better  than  did  the  average  power  detector.  For  a  false  alarm  rate  of 
4.5  false  alarms  per  hour  and  a  detection  probability  of  0.5,  the  de¬ 
tection  capability  gain  observed  was  0.17  mb  in  one  case,  and  0.25  mb  in 
the  other  case.  Shensa's  results  suggest  that  a  dominant  high-frequency 
peak  in  a  weak  explosion  signal  is  significantly  more  detectable  than  is 
the  result  obtained  when  the  power  is  averaged  over  a  broad  band  of  fre¬ 
quencies  between  0.9  and  3.6  Hz.  These  results  pertain  to  weak  explosion 
signals  which  are  nominally  4  dB  above  the  ambient  noise  level. 

As  a  step  toward  further  optimization  of  spectral  detectors, 

Shensa  suggests  using  prior  information  of  signal  and  noise  spectra  to 
perform  fixed  weighted  estimates  of  the  stacked  power  spectral  density 
for  various  observed  signal  types  (e.g.,  earthquakes,  explosions,  local 
events).  The  detector  would  then  determine  the  maximum  deflection  from 
such  weighted  power  averages  as  optimum  detectors  of  events  with  spectra 
matched  to  previously  observed  signals.  This  suggestion  is  similar  to 
that  of  the  matched  spectral  pre-filters  suggested  for  improving  the 
broadband  "Z"  detector. 

The  computer  capabilities  and  storage  requirements  needed  to  imple¬ 
ment  Shensa's  maximum  deflection  detector  are  modest,  and  can  probably  be 
accomplished  with  a  minicomputer  (distributed  processing  is  probably 
feasible  for  certain  microcomputers).  The  programming  required  is  only 
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a  slightly  expanded  version  of  the  broadband  "Z"  detector.  The  bulk  of 
that  program  handles  the  suppression  of  redundant  code  detections. 
Therefore,  the  Inclusion  of  a  set  of  front-end  filters  o\  .ront-end  FFTs 
represents  only  a  modest  program  expansion.  Since  the  maximum  deflec¬ 
tion  detector  reduces  the  spectral  data  down  to  a  single  value  of  "Z" 
at  each  point  In  time,  the  major  part  of  the  "Z"  detector  algorithm  con¬ 
cerned  with  suppression  of  redundant  coda  need  only  be  performed  once 
for  each  time  slice  of  data. 
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2.4  ANALYTIC  ENVELOPE  DETECTOR,  FARNBACH  (1975);  UNGER  (1978) 

The  mathematical  concept  behind  this  detector  is  the  representa¬ 
tion  of  the  signal  as  the  real  part  of  a  complex  function  of  time.  In 
the  simplest  case  this  amounts  to  a  phasor  representation  of  a  sinusoidal 
signal  using  a  complex  exponential  Instead  of  a  sinusoid.  In  the  case  of 
narrow  band  (but  not  purely  sinusoidal)  signals,  the  concept  Is  general¬ 
ized  to  include  amplitudes  and  frequencies  that  vary  "slowly"  with  time. 

The  advantage  of  this  representation  is  the  separation  of  amplitude  and 
phase  Information  which  in  a  real  signal  are  blended  In  a  way  which 
is  hard  to  separate  visually. 

Any  real  function  f(t)  can  be  extended  to  form  the  complex  analytical 
signal  (Figure  2.9) 

C(t)  -  f(t)  -  1  FHi(t)  (2  32) 

where 

4/  ^  (2-33) 

•  00 

is  the  Hilbert  transform.  (See,  e.g. ,  Bracewell,  1965.) 

If  we  can  assign  a  mean  frequency  w  to  the  signal,  as  we  might 
expect  for  a  narrow  v^and  seismic  signal,  then  we  can  write  the  complex 
representation  of  x(t)  as 

C  (t)eiwt  (2.34) 

where  C¥(t)  is  the  generalization  of  the  sinusoidal  phasor,  called 

A 

the  complex  envelope  function. 
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Figure  2,9.  The  complex  analytic  signal,  A  modulated  carrier  f(t), 
its  quadrature  function  F^tt),  and  the  associated  com¬ 
plex  analytic  signal  are  Si  1  shown  as  functlons’of  the 
real  variable  time  (after  Bracewell ,  1965). 


i 
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(2.35) 


Cx(t)  -  E(t)ein^ 

Then,  our  original  signal  x(t)  is  simply  the  real  part  or 
x(t)  »  E(t)  cos[co0  t  +  <fr(t)]  (2.36) 


l 


I 


where 

E(t)  Is  the  instantaneous  envelope  amplitude 
Wq  Is  the  mean  frequency 

<l>(t)  is  the  instantaneous  phase  with  respect  to  this 
mean  frequency. 

The  Instantaneous  frequency  w(t)  is  given  by 


w(t)  3  ^  <t>( t )  +  coQ 


and 


E(t)  * 


(2.37) 


(2.38) 


where  xHl- (t)  is  the  Hilbert  transform  of  x(t). 

From  a  practical  computational  viewpoint,  being  able  to  calculate 
the  envelope  function  implies  an  associated  time  window  or  epoch.  In¬ 
deed,  examination  of  Figure  2.10  shows  us  that,  in  fact,  the  envelope 
function  E(t)  (see  Eq.  (2.35))  is  simply  a  STA  of  |x(t)|.  Now, 
we  can  examine  a  phasor  plot  of  the  complex  envelope  function  for  a 
seismic  waveform  that  is  composed  of  a  signal  vector  s(t)  and  a  noise 
vector  n(t)  combining  to  form  a  resultant  phasor  x(t)  with  modulus 
|  x(t)  |  and  phase  4>(t).  If  we  assume  that  the  signal  phase  is  zero, 
the  phasor  diagrams  are  shown  in  Figure  2.11,  under  both  large  signal 
(left  frame)  and  small  signal  (right  frame)  conditions. 


» 
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Figure 


±ii  I  i 

w0  2n/w(t) 


,10,  Envelope  representation  of  the  complex  analytic  signal 
(after  Bracewell ,  1965), 
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| r 1  •  INSTANTANEOUS  AMPLITUDE 
cb  =  INSTANTANEOUS  PHASE 


Figure  2.11.  The  instantaneous  phasor  for  the  complex  analytic  signal 
(after  Unger,  1978). 
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In  practice,  of  course.  It  Is  not  possible  to  distinguish  between 
the  Instantaneous  values  of  |x(t)|  and  |n(t)|.  It  is  more  useful  to 
examine  the  probability  that  |x(t)|  >  nQ  some  long  term  property  of 
the  noise  vector  amplitude.  Unger  showed  that  If  most  of  the  time 
|n(t)|  <  nQ  then 

P(|x$(t)|  >  nQ)  <  P(|xs(t)|  >  jn(t) | ).  (2.39) 


Thus, 


P( |rs ( t) |  >  nQ)  =  e  ,  s(t)  =  0  (2.40a) 

e  <  P(|r$(t)  |  >  nQ)  <  1  -  tt'1  arccos  2[n('tir  *  0  <  1*^1 

<  |n(t) |  (2.40b) 

e  <  P(|rs(t)|  >  nQ)  <  1  ,  I s(t) |  >  2| n(t) | , 

(2.40c) 

where  the  e  is  the  probability  that  the  noise  envelope  in  the  presumed 
signal  gate  is  greater  than  nQ  in  the  lagging,  presumed  noise  gate. 
This  value  depends  on  the  statistical  distribution  of  noise  envelope 
values.  Thus,  in  the  presence  of  signal,  this  probability  is  greater 
than  e,  and  increases  with  SNR,  but  is  subject  to  an  upper  bound  which 
ranges  from  0.5  to  1.0  as  determined  by  the  instantaneous  SNR.  The  dif¬ 
ference  between  the  probability  value  and  the  upper  bound  increases  with 
the  difference  nQ  -  |n(t)|.  This  probability  distribution  function 
Is  given  in  Figure  2.12,  together  with  the  phase  bias  probability 
distribution  curve  given  later. 

Unger  (1978)  described  two  detectors,  one  based  on  the 
amplitude  of  the  envelope  function  and  the  other  based  on  its 
phase. 


r 
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PROBABILITY 


1 .  Amplitude  Detection 


Unger's  envelope  amplitude  detection  algorithm  used  simply  n^  *  jn|max 
over  a  long  term  noise  gate  and  approximated  the  probability  P(jx(t)j  >  nQ) 
by  counting  the  number  of  times  that,  in  a  presumed  signal  gate,  the  envelope 
exceeds  the  maximum  envelope  in  a  lagging  noise  gate. 

The  procedure  was  as  follows  (Figure  2.13). 

"First,  over  a  specified  warm-up  period  (e.g.,  40  seconds),  the 
peak  noise  envelope,  | n | max *  is  established.  This  peak  envelope  is  cosine 
tapered  over  subsequent  waveform  points,  with  a  spet.fied  time  constant 
(e.g.,  with  a  60-second  time  constant,  the  original  peak  value  is  halved 
at  30  seconds  and  equals  zero  at  60  seconds).  An  envelope  value  exceeding 
the  tapered  peak  value  establishes  a  new  noise  peak,  unless  a  signal 
detection  is  declared;  in  that  case  no  noise  peak  update  takes  place  until 
the  signal  is  declared  to  be  terminated. 

A  signal  detection  is  called  whenever,  in  a  forward-looking  (leading) 
time  window  of  specified  length  (e.g.,  4  seconds),  the  probability  that 
the  envelope  is  greater  than  the  tapered  peak  noise  envelope, 

P(|x  (t)|  >  |n|max)  exceeds  a  specified  threshold  TH1  (e.g.,  TH1  *  0.3). 

When  this  probability  reaches  its  maximum  the  algorithm  starts  looking 
for  the  first  signal  envelope  peak.  When  the  ratio  of  first  signal 
envelope  peak  and  tapered  noise  envelope  peak  exceeds  a  second  specified 
threshold,  the  SNR  threshold  TH2  (e.g.,  TH2  =  2  to  3  dB) ,  the  signal 
detection  is  confirmed  and  a  frequency-dependent  stepback  is  performed 
to  determine  the  signal  onset  time. 

The  stepback  procedure  (Figure  2.14)  is  based  on  the  observation 
that  in  most  cases  the  first  signal  envelope  peak  (at  t^)  occurs  within 
one  signal  period,  and  frequently  at  approximately  3/4  period,  after  the 
signal  onset  (at  tQ).  In  a  high-SNR  waveform  the  signal  onset  time  is 
most  accurately  found  by  detecting  the  first  maximum  or  minimum  of  the 
signal's  instantaneous  value  (at  t^),  and  stepping  back  1/4  period 
(=  0.25/instantaneous  frequency  at  t3).  For  low-SNR  waveforms  the  first 
quarter  period  may  be  obscured  by  noise;  in  that  case  we  step  back  3/4 
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mean  period  (=  0.75/mean  frequency  at  t^)  from  the  first  signal  envelope 
peak  at  t^.  The  mean  frequency  is  the  closed-form  derivative  of  the 
phase  regression  polynomial  evaluated  at  time  t^.  The  search  for  the 
first  quarter  period  is  started  at  t2,  i.e.,  at  0.8  mean  period  before 
t^;  the  first  quarter  period  is  detected  when  its  maximum  or  minimum 
exceeds,  by  a  third  threshold,  TH3  (e.g.,  TH3  =  1  dB),  the  immediately 
preceding  noise  in  the  one-second  time  interval  (tj.  t?). 

If  the  second  threshold  (the  SNR  threshold)  is  not  satisfied,  the 
detection  is  annulled  and  the  noise  peak  value  is  updated  with  what  at 
first  was  believed  to  be  the  signal  envelope  peak.  Thereafter,  the  noise 
peak  is  updated  as  usual,  until  the  next  supposed  signal  detection,  etc. 

The  signal  end  time  is  found  as  the  moment  of  the  first  envelope 
minimum  occurring  either  after  P(|rs(t)|  >  |n|  )  falls  below  its 

threshold,  or  after  the  signal  duration  exceeds  a  specified  maximum, 
whichever  is  first.  If  this  envelope  minimum  is  greater  than  the  tapered 
noise  peak  the  noise  peak  envelope  is  updated  with  this  value,  and  noise 
peak  updating  and  signal  detection  resume  as  normal.  This  procedure 
enables  the  detection  and  timing  of  later  phases  and  other  signals  in 
the  coda." 

Sax  (1980)  reported  that,  "in  its  present  state,  Unger's  detector  op¬ 
erated  on  a  single  trace  at  a  70  percent  detection  probability  with  about 
seven  false  alarms  per  hour.  Its  detection  performance  is  certainly  no 
greater,  and  is  probably  less  than,  that  of  the  other  detectors  discussed 
here.  Nonetheless,  it  has  demonstrated  superior  performance  in  accurately 
timing  P-wave  signals. 

Unger  (1978)  tested  his  detector  against  an  analyst's  timing  of 
P-wave  signals.  Nearly  half  of  the  events  examined  were  timed  with  no 
apparent  error.  Ninety  percent  of  the  signals  were  timed  within  +  0.5 
seconds  of  the  analyst  pick.  Slowly  emergent  earthquake  signals  were 
sometimes  picked  several  seconds  late.  The  standard  deviation  for 
arrival  times  determined  using  Unger's  detector  is  0.2  seconds.  By 
comparison,  the  "Z"  Detector  times  signals  1.70  seconds  late  on  an 
average,  with  a  standard  of  2.20  seconds.  Thus,  Unger's  detector  demon¬ 
strated  superior  capabilities  for  automatically  timing  seismic  signals. 
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As  a  result  of  Its  superior  capability  in  timing  seismic  .-'  waves, 
Unger's  detector  was  used  in  the  1979  VSC  Event  Identification  Experiment 
by  ENSCO  to  automatically  time  and  edit  P-wave  signals  from  long  seismic 
records.  For  example,  a  five  minute  record  containing  a  possible  P-vave 
signal  is  automatically  processed  by  Unger's  detector  to  transform  it 
into  a  50-second,  signal entered  edit  of  the  P  wave. 

The  experience  with  this  application  was  that  almost  all  of  the 
detected  sicrals  were  accurately  timed,  with  only  a  negligible  number  of 
false  alarms.  However,  problems  were  encountered  with  missed  signals. 

These  included  both  impulsive  signals  of  very  short  duration  and  gradually 
emergent  signals.  The  missed-signal  problem  could  be  corrected  by 
employing  variable-length  time  gates  for  the  forward-looking  signal  window 
and  by  using  ordered  noise  statistics  rather  than  the  maximum  noise 
estimates  employed  by  Unger.  In  a  few  cases,  Unger's  detector  inadvertently 
shut  itself  off  when  it  encountered  large  glitches  or  spikes  in  the  noise 
preceding  a  signal.  Those  cases  could  have  been  avoided  by  using  robust 
median  estimates  of  noise  rather  than  the  maximum  noise  estimates  used  by 
Unger. 

2.  Phase  Detection 

Unger  also  investigated  the  probability  distribution  for  the  phase 
<j>(t).  For  a  noise  phase  angle  uniformly  distributed  between  ±tt,  in  the 
absence  of  a  signal  P ( | d> ( t)  |  <  rr/2) ,  equals  0.5.  "The  vector  diagrams  in 
Figure  2.11  show  that,  when  a  signal  is  present,  the  phase  angle  is 
statistically  biased,  i.e.,  the  above  probability  is  greater  than  0.5. 

For  a  given  instantaneous  SNR,  this  probability  is: 


’(lt(t)|  <  ir/2)  =  0.5 


,  s(t)  =  0 


(2.41a) 


P( 1 4>(t) |  <  ir/2)  =  1  -  it"1  arccos  liHH  ,  |s(t)|  5  n(t) 

I n( t) | 


(2.41b) 


P(U(t)!  <  it/2)  =  1 


,  I  s( t)  I  >  I  r»( t)  I  .  (2.41c) 


This  is  plotted  in  Figure  2.12. 
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For  a  constant  SNR,  this  phase  bias  probability  can  be  approximated 
by  counting,  within  a  waveform  window  of  sufficient  length,  the  number  of 
times  that  the  phase  fluctuation  Is  within  +  tr/2  radians,  and  dividing  this 
count  by  the  number  of  window  points.  However,  since  In  general  the  SNR 
will  vary  Inside  the  window,  this  approximated  phase  bias  probability  will 
not  relate  to  the  SNR  exactly  as  in  Eqs.  (2.44).  Nevertheless,  it 
still  is  some  measure  of  the  average  SNR  in  the  window,  and  is  a  good 
detection  parameter,  since  it  will  have  a  value  greater  than  0.5  if  a 
signal  is  present  in  the  window.  In  this  manner,  phase  detection  is 
established  in  principle.  In  the  application  to  actual  data,  however, 
there  are  some  complications  which  will  be  discussed  shortly. 

The  phase  bias  probability  distribution  function  is  compared  to 
the  envelope  detection  probability  distribution  function  in  Figure  2.12. 

We  observe  the  important  fact  that  the  phase  distribution  curve  reflects 
a  detection  sensitivity  which  is  twice  that  of  the  envelope,  since  the 
arccosine  argument  equals  the  instantaneous  SNR  in  the  case  of  phase 
detection,  and  only  one  half  the  instantaneous  SNR  in  the  case  of  envelope 
detection.  This  suggests  that,  in  principle,  phase  detection  is  at  least 
6  dB  better  than  envelope  detection,  especially  when  regarding  the  fact 
that  the  envelope  curve  represents  the  upper  bound  of  envelope  detection 
sensitivity.  The  detection  sensitivity  of  the  instantaneous  phase  has 
been  shown  and  used,  for  instance,  in  underwater  sound  propagation  studies 
(Steinberg  and  Birdsall  (1966);  Unger  and  Veenkant  (1967a,  1967b)."  In 
the  work  of  Unger  (1978)  this  technique  was  applied  to  the  detection  and 
timing  of  seismic  signals. 

Since,  in  general,  the  signal  phase  varies  with  time  in  a  determinis¬ 
tic  manner  (e.g.,  in  LP  dispersed  waveforms),  the  principle  of  phase  detec¬ 
tion  can  only  be  applied  in  those  cases  where  a  model  for  the  expected 
signal  phase  angle  variations  can  be  adequately  specified.  Clearly,  such 
a  model  can  be  specified  for  most  LP  waveforms,  but  it  is  not  so  obvious 
for  most  SP  waveforms.  Also,  cor.  ary  to  the  assumption  above,  the  noise 
phase  is  not  uniformly  distributed,  but  may  rather  follow  a  somewhat 
more  deterministic  trend.  This  is  the  case,  for  instance,  when  the 
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dominant  noise  frequency  differs  from  the  reference  frequency,  f q,  thus 
causing  a  linear  phase  trend.  The  "continuous"  phase  then  may  traverse  a 
number  of  cycles  within  a  given  time  gate.  These  facts  necessitate 
"tracking"  the  Instantaneous  phase  function;  the  phase  fluctuations  about 
the  tracked  or  time-variant  mean  phase  then  may  be  studied  for  signal 
detection.  For  noise,  the  fluctuations  should  be  randomly  distributed;  in 
the  presence  of  signal  they  will  be  statistically  biased.  Thus,  the 
performance  of  the  phase  detector  now  rests  with  the  efficiency  of  the 
tracking  process  with  respect  to  some  presumed  model  governing  the  phase 
variations  of  signals,  and  also  with  the  validity  of  the  model  used  to 
estimate  the  signal  phase  angle." 

Unfortunately,  the  rapid  change  of  phase  associated  with  a  body 
wave  arrival  means  that  few  sample  points  are  available  for,  say,  regression 
analysis  of  the  phase  changes,  unless  the  data  is  very  oversampled.  This 
inability  to  track  phase  led  Unger  to  abandon  phase  detection,  the  theoreti¬ 
cally  more  sensitive  detector,  and  rely  simply  on  envelope  detection. 

2.5  ALLEN  DETECTOR,  ALLEN  (1978) 

This  detector  is  based  upon  a  heu.’istic  detection  statistic 

E(t)  -  [x(t)]2  +  C  (2.45) 

where  x(t)  is  the  input  signal.  The  implementation  of  this  proceeds 
as  follows: 

1.  Calculate 

AYi  =  xi  -  Vi  •  (2.46) 

2.  Calculate 

Yi  =  aYi-i  +  AYi  •  (2.47) 

3.  Calculate  the  "statistic" 

E.  =  +  C  AY*  .  (2.48) 
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We  can  think  of  the  detection  statistic  as  being  formed  from 
the  weighted  sum  of  the  "power"  of  two  filter  outputs,  one  being  a  high- 
pass  (really  a  bandpass)  and  the  other  a  differentiator,  as  shown  in 
Figure  2.15.  The  second  filter  clearly  accentuates  the  energy  near  the 
Nyqulst  frequency,  making  good  arrival  time  estimates  possible. 

The  STA  and  LTA  of  the  detection  statistic  E  are  calculated  in 
the  normal  way  with  recursive  filters: 

STA1  =  (1  -  C$)STA1_1  +  C$  E.  (2.49) 

LTA.  =  (1  -  Cl)LTAm  +  CL  E.  (2.50) 

Finally,  an  event  is  declared  if  the  STA/LTA  ratio  exceeds  some  thres¬ 
hold  value. 

The  rest  of  Allen's  algorithm  is  concerned  with  verifying  that  the 
declared  event  is  not  a  false  alarm  and  with  timing  the  arrivals.  The 
program  p»*oceeds  as  follows: 

1.  When  an  event  is  declared  the  time  is  recorded  along 
with  Y  and  the  first  difference  of  Y. 

2.  "The  program  now  enters  a  pair  of  nested  loops  in  which  it 
searches  for  a  peak  amplitude  and  the  subsequent  zero 
crossing.  When  the  zero  crossing  is  detected,  the  program 
breaks  out  of  the  loop,  records  the  zero-crossing  time 
relative  to  event  onset  and  the  signed  amplitude  of  the 
preceding  peak.  The  time  and  amplitude  information  is 
stored  ...  for  later  use  by  analysis  routines"  (Figure  2.16). 

3.  At  each  zero  crossing  of  Y,  the  value  of  the  STA  is  compared 
with  a  constant  and  a  record  kept  of  the  number  of  timts,  S, 
that  the  STA  is  less  than  this  constant. 

4.  The  "event  will  be  declared  over  when  some  number  of  zero 
crossings  with  the  STA  less  than  the  constant"  have  occurred. 

This  termination  number,  L,  will  be  small,  typically  3,  at 
the  start  of  an  event  to  enable  the  program  quickly  to  reject 
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Figure  2.15.  (a)  Block  diagram 

sponse  functions  < 


the  Allen  detector  and  (b)  the  re- 
the  two  filter  elements. 


A* 


» 


Figure  2.16.  Schematic  earthquake  with  data  stored  during  events. 

(a)  Seismic  event  with  onset  and  end  points,  (b)  Dots 
represent  zero-crossing  times  and  previous  peak  ampli¬ 
tudes.  Amplitude  bars  indicate  background  noise  preced¬ 
ing  onset,  and  arrow  gives  first  difference  at  onset. 

(c)  Data  stored  for  use  by  analysis  routines  after  event 
terminates  (after  Allen,  1978). 
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noise  spikes  or  other  very  short-term  Interference.  When 
the  algorithm  is  well  into  a  larger  event,  however,  l  must 
be  considerably  larger  to  ensure  chat  an  earthquake 
observation  is  not  terminated  too  early  during  a  quiet 
period  between  phase  arrivals.  The  present  version  of  the 
program  uses  the  relation 

L  -  3  +  M/3  (2.51) 

where  M  Is  the  current  number  of  observed  peaks  in 
an  event. 

5.  This  decision  of  whether  the  event  is  ove^  is  simply  a 
comparison  of  S  and  l,  with  branches  to  continue  the 
program  in  the  observation  loop  or  tc  terminate  observa¬ 
tion  of  the  event  as  required. " 

6.  A  test  is  made  for  the  duration  of  the  event  to  eliminate 
such  noise  bursts  as  line  spikes,  etc.  For  small  local 
events,  Allen  used  typically  the  criterion  that  the  event 
should  be  longer  than  1.5  seconds  and  have  recorded  more 
than  40  peaks. 

The  performance  statistics  currently  available  on  Allen's  detector 
are  based  on  local  seismic  network  data  sampled  200  times  a  second.  The 
results  obtained  are  extrapolated  below  to  those  for  a  20-Hz-sampled  regional 
or  teleseismic  network. 

The  operating  characteristics  demonstrated  by  local  network  opera¬ 
tion  show  a  70  percent  detection  of  analysts'  picks,  with  36  false  alarms 
over  44  hours  of  operation.  Allen's  detector  also  grades  the  quality  of 
signal  detections.  None  of  the  36  false  alarms  as  graded  at  the  highest 
level  of  reliable  detections. 

Allen's  detector  times  local  event  signals  to  a  standard  deviation 
of  0.05  seconds  with  200-Hz  local  network  data.  This  would  scale  to 
0.5  seconds  for  our  20-Hz  data,  much  better  than  the  2.2  second  standard 
deviation  of  the  "Z"  Detector. 


,<.>.flia’ima a  rrr. 
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In  comparison  with  the  "Z"  Detector  at  the  level  of  70  percent 
detection  probaoility,  the  Allen  Detector  operated  with  0.82  false  alarms 
per  hour  as  a  local  event  detector  compared  with  0.4  false  alarms  per  hour 
of  DET  459  and  with  0.8  false  alarms  per  hour  at  KSRS.  By  design,  the  "Z" 
Detector  is  expected  to  maintain  a  more  stable  level  of  false  alarm  control. 
Allen's  detector  concept  applied  to  the  detection  of  regional  and  tele- 
seismic  signals  sampled  at  20  Hz  is  now  being  studied  by  ENSCO's  DSC 
Division. 

It  is  interesting  that  Allen's  detector  operates  on  a  completely 
different  principal  than  does  the  "Z"  Detector.  The  power  of  Allen's 
detector  stems  from  the  sophisticated  post-detection  analysis  of  possible 
signals  which  is  designed  to  eliminate  most  false  alarms.  The  frequency 
is  estimated  by  counting  zero  crossings.  Also,  the  time  duration  of  the 
signal  is  estimated.  A  minimum  event  duration  of  1.5  seconds  and  an 
acceptable  range  of  frequencies  are  some  post-detection  analysis  require¬ 
ments  for  accepting  detections  as  representing  valid  signals. 

There  appears  to  be  room  for  a  substantial  performance  gain  •'n 
detecting  P  waves  by  using  the  Allen  Detector,  and  especially  by  optimizing 
the  post-detection  analysis  for  that  purpose.  The  detector  requires  a 
limited  memory  (about  1.0  K  bytes)  and  a  very  modest  computer  capability. 

The  characteristic  function  and  post-detection  processing  of  the  Allen 
Detector  is  in  a  form  most  suitable  for  optimizing  and  training  the 
automatic  detector  to  track  the  performance  of  experienced  analysts. 
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2.6  THE  STEWART  TETECTOR,  STEWART  (1977) 


This  detector  uses  an  algorithm  that  was  designed  to  mimic  the 
mental  process  tnat  a  trained  analyst  performed  on  a  seismic  trace.  A 
nonlinear  high  pass  filter  is  applied  to  the  signal  to  transform  it  in 
such  a  manner  that: 

a.  The  oscillatory  nature  of  the  signal  is  preserved; 

b.  Th'  direction  of  the  first  motion  is  preserved; 

c.  It  high  passes  the  signal  in  the  normal  manner,  thus  re¬ 
ducing  the  D.C.  and  drift  components; 

d.  Slightly  emergent  onsets  of  seismic  signals  are  enhanced. 

The  transformation  process  is  illustrated  in  Figures  2.17  and  2.18. 
The  "algorithm  computes  the  modified  seismic  signal  MDXk  from  the  in¬ 
coming  signal  Xk,  where  k  represents  the  current  time  epoch.  The  re¬ 
sults  of  applying  this  transformation  to  a  1-Hz  and  a  3-Hz  sinusoidal  signal 
are  shown  in  Figure  2.18.  From  this  figure  the  high-pass  nature  of  the 
transformation  is  apparent. 

The  first  step  in  the  transformation  process  (Figure  2.17)  is  to 
compute  the  simple  first  difference  of  the  incoming  signal  (DXk  =  Xk  - 
Xk_.|  where  k  represents  the  current  time  epoch).  The  sign  of  the 
current  first  difference  DXk  is  compared  to  the  sign  of  the  previous 
first  difference  DXk_-|.  If  the  signs  are  the  same,  and  if  this  sign  has 
persisted  for  less  than  eight  consecutive  times,  then  the  value  of  the 
modified  signal  MDXk  is  taken  to  be  its  current  value  increased  by  DXk> 
Otherwise,  the  value  of  the  modified  signal  is  taken  to  be  DXk."  This 
technique  transforms  the  incoming  signal  in  such  a  way  that  the  four  ob¬ 
jectives  listed  above  are  accomplished.  "This  modified  seismic  wave 
form  is  the  basis  for  nearly  all  signal  analysis  by  the  on-line  system. 

T^e  one  exception  to  this  is  that  the  maximum  signal  amplitude  is  determined 
from  the  incoming  seismic  signal  Xk.!l 


j'AtwirfnMiiHw, 
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Figure  2.17.  Processing  steps  to  compute  the  modified  seismic  signal 

MDX*  from  the  incoming  signal  X^.  The  modified  signal  is 
used  extensively  in  the  detection  logic  of  the  on-line 
system  (after  Stewart,  1977). 
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Figure  2.18.  Effect  of  the  filtering  operation  (summarized  in  Figure 
2.17)  on  a  1  Hz  and  3  Hz  sinusoidal  signal.  The  upper 
signal  is  the  unfiltered  analog  output  of  a  function 
generator,  fed  simultaneously  to  one  channel  of  the  on¬ 
line  system  and  to  a  chart  recorder.  The  lower  signal 
is  the  filtered  function  MDX|<  computed  in  real-time  by 
the  on-line  system  and  routed  through  a  digital-to- 
analog  converter  to  the  chart  recorder.  Digitizing  rate 
is  50  samples/second.  Peak  amplitude  of  the  3  Hz  filtered 
signal  is  approximately  twice  that  of  the  1  Hz  signal 
(after  Stewart,  1977). 
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The  long-term  averages  of  three  quantities  are  computed; 

1.  The  mean  of  the  signal  itself 

LTA1.  =  (1  -  Cl)  LTAli_1  +  Cl  X.  . 

2.  An  estimate  of  the  standard  deviation  of  the  signal  is  made 
by  calculating  the  LTA  of  the  absolute  value  of  the  varia¬ 
tions  of  the  signal  about  its  mean 

LTA2i  =  (1  -  C2)  LTA2._^  +  C2  |X.  -  LTA1.|  . 

3.  A  similar  estimate  of  the  standard  deviation  of  the  modified 
signal  is  the  LTA  of  the  absolute  value  of  the  MDX^ 

LTA3.  *  (1  -  C3)  '-TA3i_1  +  C3  |MDX.| 

where  Cl,  C2  and  C3  are  the  filter  constants. 

These  quantities  are  intended  to  characterize  the  signal  under  the 
no-event  hypothesis  and  thus  updating  their  values  is  suspended  if  a 
tentative  event  is  detected.  The  detection  algorithm  is  a  two- stage  pro¬ 
cess  with  the  first  stage  detection  test  being  simply  a  test  of  whether 
MDX^./LTAS^  is  larger  than  some  constant  C4. 

When  a  tentative  detection  is  made,  the  algorithm  proceeds  to  the 
post-detection  stage  designed  to  reduce  the  FAR.  Stewart  (1977)  divides 
this  process  into  two  modes,  the  P-phase  processing  and  the  coda  pro¬ 
cessing. 

The  P-phase  processing  continues  for  0.5  second  after  the  tentative 
detection  is  declared.  (The  data  he  considered  were  sampled  at  50  sps.) 
During  this  interval,  the  following  four  criteria  must  be  met  to  proceed 
or  else  the  tentative  detection  is  cancelled: 


1. 


MDX1 

LTA3T 


>  C4 


N1  times; 
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2.  The  time  duration  for  which 
MDX. 

rjj -1-  <  C4  must  not  exceed  0.28 

L  MJi  seconds  consecutively; 

3.  (max  j  >  C5 
max  X .  j 

4’  L?A2”  >  8  * 

If  these  four  tests  are  passed,  the  algorithm  then  proceeds 
with  the  coda  processing.  This  process  continues  for  a  time  interval 
determined  by  testing 

>  C6  , 

if  this  test  fails  for  a  continuous  2  seconds,  the  coda  processing  is 
terminated.  The  following  three  tests  must  be  passed  during  this  phase 
of  the  algorithm. 

1.  Coda  length  must  be  _>  4  seconds; 

MDX.  MDX. 

2'  TJA3T  >  06  and  LTA3T  <  ‘  C6 

*  J 

is  alternating  sequence  at  least  six  times; 

3.  The  number  of  oscillations  of  MDX..  must  exceed  0.5  H2/time 
duration  of  coda. 

This  algorithm  was  implemented  by  the  USGS  using  data  from  their 
central  California  and  Oroville  networks.  Stewart  (1977)  reports  on  test 
results  from  this  network  obtained  during  a  one  month  period  which  in¬ 
cluded  a  large  aftershock  sequence  near  Oroville. 

Table  2.1  (after  Stewart,  1977)  summarizes  the  results  of  both 
detection  and  on  line  location  for  the  Oroville  network  while  Table  2.2 
(also  from  Stewart,  1977)  summarizes  the  results  for  one  month  of  opera¬ 
tion  on  the  central  California  network. 
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TABLE  2.1 


DETECTION  AND  LOCATION  CAPABILITY  OF  REAL-TIME  SYSTEM  FOR  THE 
OROVILLE  SEISMIC  NETWORK  DURING  OCTOBER  1975 


of  Events 

Percentage 

Comment 

107 

90.7 

Detected  and  located. 

8 

6.8 

Detected  but  not  located,  because 
data  were  too  noisy  or  too  sparse. 

1 

0.8 

Not  detected. 

2 

1.7 

Not  detected,  preceded  within  60 
seconds  by  another  e/e nt 

118 

100.0 

I 


4 


TABLE  2.2 

DETECTION  CAPABILITY  OF  REAL  TIME  SYSTEM  FOR  THE  CENTRAL 
CALIFORNIA  SEISMIC  NETWORK  DURING  OCTOBER  1975 


No.  of  Events  Percentage 


225 

86.5 

Detected. 

25 

9.6 

Not  detected,  events  north  of  38°30'. 

5 

1.9 

Not  detected. 

3 

CM 

rH 

Not  detected,  computer  maintenance  in 
progress. 

2 

0.8 

Not  detected,  preceded  within  60  seconds 

260 

100.0 

by  another  event. 
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A  portable  field  recorder.  Teledyne  Geotech's  Model  MCR-3G0 
Microcorder  {Veith,  1978),  which  uses  this  algorithm,  has  been  developed 
and  tested.  The  tests  were  run  on  two  field  units  operating  at  the 
Nevada  Test  Site  where  the  noise  sources,  both  natural  and  cultural, 
varied  both  in  character  and  time.  Oata  was  sampled  at  rates  varying  from 
20  sps  to  200  sps  for  several  days.  Comparison  between  the  results  of 
this  algorithm  and  the  "standard  STA/LTA"  algorithm  was  made.  The  main 
results  of  these  tests  were: 

1.  "Once  triggered,  the  STA/LTA  detector  often  shut  down  during 
a  relatively  low  signal  period  only  to  trigger  again  after 

a  second  or  two.  The  algorithm  detector  did  not  shut  down 
until  the  end  of  the  signal. 

2.  An  "emersio"  signal  would  not  trigger  the  algorithm  unless  it 
contained  a  higher  amplitude  energy  burst.  The  STA/LTA  de¬ 
tector  often  triggered  from  "emersio"  signals.  (Emersio  is 
used  here  in  the  sense  that  both  the  long  and  the  short  term 
averages  are  increasing  slowly  in  the  STA/LTA  ratio  is  near 
the  specified  SN  value.)" 

Veith  (1978)  finally  concludes  that,  "tests  on  limited  data  show 
the  algorithm  can  ignore  noise  while  detecting  earthquake  energy.  It 
can  make  the  transition  between  changes  in  noise  level  with  every  little 
difficulty  and  record  seismic  signals  within  virtually  any  background 
without  changing  the  settings  of  the  algorithm  parameters.  The  algorithm 
will  record  complex  signals  without  the  termination  problems  typical 
of  the  STA/LTA  detector. 

While  the  tests  gave  very  good  results  in  detecting  seismic  events, 
they  were  not  perfect.  It  may  be  possible  to  develop  additional  simple 
tests  during  the  field  test  period  which  would  provide  a  more  powerful 
discriminant.  The  principle  difficulty  lies  in  obtaining  tests  which 
may  be  passed  by  microearthquake  signals,  teleseisms  or  surface  waves  in 
accordance  with  the  pass  band  selected." 
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2.7  WALSH  DETECTOR,  GOFORTH,  T.  AND  E.  HERRIN  (1980) 


This  detector  is  based  on  the  representation  of  the  digital  seismic 
signal  In  terms  of  Walsh  functions.  These  form  an  ordered  set  of  rec¬ 
tangular  waveforms  which  are  either  +1  or  -1.  They  are  ordered  ac¬ 
cording  to  the  number  of  zero  crossings  per  interval.  Like  the  exponential 
functions  of  the  Fourier  transform,  they  form  an  orthonormal  set  over 
some  interval.  Thus,  the  sampled  seismic  signal  can  be  represented 
in  terms  of  its  Walsh  transform  coefficients  by, 


N-l 

xj  =2Wk  Wa1(k’j)  (2.52) 

k*0 

where 

N-l 

Wk  3  1^2  Xj  Wal(k.j)  (2.53) 

j=0 


and  k  is  termed  the  sequency.  Figure  2.19  illustrates  the  Walsh  functions 
Wal(N.i)  for  N  *  0,  ...»  8.  This  transformation  may  be  compared  to  the 
discrete  Fourier  transform  pair  of: 


x 


j 


eiwkj 


(2.54) 


where 


N-l 

E 

j=o 


xje 


-iojjk 


and 


2tt 

N 
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Figure  2.19.  The  first  eight  Walsh  functions. 
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Obviously,  since  Walsh  transforms  are  generated  using  simple 
rectangular  +  or  -  unity  functions  rather  than  sines  and  cosines, 
they  may  be  calculated  much  more  efficiently  (quicker)  in  a  computer. 
Clearly,  to  generate  the  Walsh  transform  of  the  discrete  x^  all  that 
must  be  done  are  logical  compares  and  adds.  Performing  functions  such  as 
prewhitening  and  filtering  will  similarly  be  much  more  efficient  than 
their  time  domain  or  Fourier  transform  counterparts.  Indeed,  the  purpose 
of  this  detector  is  speed  or  computational  efficiency  in  the  implementa¬ 
tion  of  concepts  used  In  other  detectors. 

A  block  diagram  of  the  detection  algorithms  is  shown  in  Figure 
2.20.  The  digital  data  is  analyzed  in  64  sample  windows  (3.2  seconds  for 
20  sps  data)  with  a  32  sample  overlap.  The  analysis  proceeds  as  follows: 


1.  The  Walsh  transform  of  the  64  data  samples  is 
calculated 


63 

Wk(x)  =  ^:^Xj  Wa1(k’j)  k  *  °*  ••••  63  (2.57) 

j=0 


2.  The  Walsh  coefficients  are  multiplied  by  prewhitening  weights 
Nr  calculated  beforehand  (see  below). 

3.  Weights  Bk  of  0  or  1  are  applied  to  isolate  the  fixed 
signal  band. 

4.  The  absolute  values  of  the  weighted  Walsh  coefficients 
are  summed  to  give  the  final  detection  "statistic" 


8k-1  klikik2 
=  0  for  other  k 


(2.58) 

(2.58) 


The  detection  statistic  Y  is  thus  equivalent  to  the  STA  here  calculated 
over  at  3.2  second  gate  for  20  sps  data. 
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The  prewhitening  weights  Nk  are  currently  calculated  In  a  non- 
adaptive  manner  as  follows: 

1.  A  34  minute  sample  of  noise  is  selected. 

2.  The  Walsh  transform  of  640  consecutive  64  sample  windows 
is  calculated. 

3.  The  mean  of  the  absolute  values  of  each  of  the  64  Walsh 
coefficients  Is  computed. 

?  4.  Weights  of  the  form  2  •  n,  where  n  is  an  integer, 

are  selected  to  whiten  the  means  in  each  subsequency 
band. 

The  detection  threshold  against  which  the  statistic  Y  is  tested 

*  is  calculated  from  the  median  value  Y50’  and  the  75  percentile  value 
Y^5  of  the  previous  512  values  of  the  STA  (14  minutes)  by  the  weighted 
combination. 

•  Threshold  =  (1  -  k'  Y5Q  +  k  Y?5  .  (2.59) 

The  coefficient,  k,  is  typically  five.  Note  that  this  is  not  the  same 

as  the  action  of  a  recursive  filter  on  Y.  If  Y  is  Gaussian  distributed 

2 

I  with  mean  y  and  variance  a  then 

Y50  =  yY 


Y75  =  1,5  aY  * 


Thus,  in  this  case. 


Threshold  =  yy  +  5  (1.15  ay  -  Uy)  . 


(2.60) 


If  the  current  value  of  the  sum  of  the  absolute  values  of  the  Walsh 
coefficients  exceeds  the  threshold  for  two  consecutive  time  windows, 
an  event  is  declared.  If  it  does  not  exceed  the  threshold,  the  sum  of 
the  absolute  values,  Y,  is  ranked  among  the  previous  512  values,  the 
oldest  value  being  discarded.  In  this  way  an  adaptive  detection  threshold 
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Is  maintained,  the  adaptation  window  being  approximately  14  minutes. 

If  a  signal  is  called,  the  threshold  is  not  updated. 

Three  separate  tests  of  the  Halsh  detector  have  been  reported  by 
Goforth  and  Herrin  (1980).  Two  consisted  of  tests  conducted  on  simulated 
data  constructed  of  no.,e  plus  scaled  8  second  long  segments  of  a 
Novaya  Zemlya  explosion  added  at  specific  t'lmes  to  the  noise  samples. 

The  noise  samples  were  taken  from: 

1.  The  SMU  KS36000  seismometer  at  Dallas  with  the  signals 
added  at  a  SNR  of  1.0  and  0.75. 

2.  The  ANMO  SRO  with  the  signals  added  to  scale  a  4.5  m^ 
and  4.2  m^  event. 

The  results  of  these  tests  were: 

la.  80  detections  (100%)  0  false  alarms 

b.  40  detections  (100%)  1  false  alarm 

2a,  39  detections  (97.5%)  0  false  alarms 

b.  36  detections  (90%)  1  false  alarm 

Finally,  thr  detection  algorithm  was  run  on  the  center  element, 
of  one  of  the  NORSAR  subarrays  for  a  five  hour  period  of  data.  The  re¬ 
sults  were  a  seven  out  of  eight  detection  success  (missing  a  3.6  m^  at 
6,222  km  epicentre!  distance)  with  a  FAR  of  0.6/hour.  In  another  over¬ 
lapping  seven  hour  test,  the  algorithm  succeeded  in  detecting  nine  out 
of  ten  events  (including  a  local  teleseism  that  the  NORSAR  beams  missed) 
with  a  FAR  of  0.86. 
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2.8  MARS  DETECTOR,  MASSO,  et  al_.  (1979) 

The  MARS  (Multiple  Arrival  Recognition  System)  detector  uses 
quasi -harmonic  decomposition  to  analyze  a  broadband  signal  by  passing 
it  through  many  narrow  band  filters.  One  thus  obtains  a  time- frequency 
breakdown  of  the  non-stationary  signal  power.  An  "event"  consists  of 
a  statistically  significant  fluctuation  from  the  random  background 
pattern  of  the  instantaneous  spectrum.  There  are  certain  strong 
affinities  between  the  MARS  detector  and  the  deflection  detector  de¬ 
scribed  in  Section  2.3,  particularly  the  "average  deflection"  (3)  ver¬ 
sion.  Both  use  frequency  domain  methods,  and  both  search  for  short 
term  increases  in  power  over  an  adaptive  frequency  band.  Unlike  the 
deflection  detector,  however,  the  MARS  detector  dees  not  calculate 
directly  a  statistic,  but  rather  contains  a  highly  nonlinear  and  adap¬ 
tive  pattern  recognition  algorithm  which  seeks  undispersed  alignments  of 
peaks  in  the  narrow  band  envelope  functions.  The  function  of  the  STA,  the 
fundamental  time  Tq,  is  roughly  dependent  on  the  Q  of  the  filters,  for 
the  narrower  the  frequency  resolution,  the  larger  the  time  duration  of 
the  impulse  response.  There  is  also  inherent  in  the  method  an  LTA,  for 
envelope  maxima  only  have  meaning  to  the  extent  they  deviate  from  past 
behavior.  But  even  though  these  two  concepts  carry  through  in  the  MARS 
detector,  they  do  not  appear  in  a  mathematically  tractible  form,  since 
the  final  decision  is  based  principally  on  a  band  width  and  dispersion 
criterion. 

The  concept  of  a  detector  based  upon  multiple  bar  J  pass  filters 
goes  back  at  least  to  Moltshan,  et  al_.  (1964),  who  stated: 

"The  signal  from  the  seismogram  is  passed  through  a  group 
of  linear  filters.  All  filters  are  divided  into  several, 
in  our  case,  two,  sub-groups.  In  each  sub-group  the  filters 
have  the  same  pass  band,  but  difference  resonance  frequencies; 
their  pass  bands  are  not  overlapping  and  cover  the  whole 
range  of  frequencies.  Then  the  output  of  each  filter  is  com¬ 
pared  with  some  threshold  which  is  chosen  according  to  the 
level  of  noise  in  this  filter.  If  in  at  least  one  of  the 
filters  the  threshold  is  surpassed,  then  the  presence  of 
a  useful  signal  —  tue  alarm  —  is  announced.  The  record 
may  begin  slightly  before  the  alarm;  which  is  why  a  delay 
line  is  desirable." 
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The  potential  advantages  in  this  technique  lie  in  the  area  of 
adaptibility  and  in  signal  characterization.  This  latter  aspect,  indeed, 
was  powerfully  utilized  in  the  work  of  Masso,  et  aH.  (1979)  which  in¬ 
volved  both  detection  and  discrimination  of  seismic  signals.  They 
describe  the  detection  part  of  the  MARS  processor  as  follows: 

"The  central  feature  of  the  detector  module  is  the  use  of  a  set  of 
narrow-band  frequency  filters  to  break  up  or  decompose  a  time  series  con¬ 
sisting  of  signal  plus  noise  into  a  set  of  ^uasi -harmonic  modulated 
"signals."  This  set  of  filtered  signals,  one  for  each  filter  of  center 
frequency  f  ,  can  then  be  used  to  determine  the  energy  arrival  time 

L 

(or  group  arrival  time,  tg)  and  amplitude  of  the  original  broad-band 
signal  by  analysis  of  the  time  modulation  of  the  filter  outputs. 

Further,  both  the  instantaneous  phase  and  frequency  of  the  individual 
filter  outputs,  that  is  the  apparent  phase  and  frequency  of  the  quasi - 
harmonic  filter  outputs  as  a  function  of  time,  can  be  determined  quite 
simply.  Thus,  the  decomposition  of  the  original  signal  wavetrain, 
possibly  composed  of  many  individual  signal  pulses,  into  quasi-harmonic 
signals  provides  the  means  of  determining  arrival  time,  amplitude  and 
phase,  all  as  functions  of  frequency.  This  then  is  the  basic  signal 
information  that  can  be  used  to  detect  a  given  type  of  signal  in  terms 
of  its  dispersion  characteristics  and  to  obtain  its  spectrum  as  well 
as  its  time  and  amplitude  relationship  with  respect  to  other  signals 
present  in  a  complex  wavetrain. 

The  basic  approach  is  to  identify  patterns  from  the  signal  and 
noise  information  as  it  is  expressed  in  the  t  -  t  plane  (the  group 
arrival  time,  tg,  versus  frequency  plane).  The  pattern  to  be  searched 
for  in  the  t  -  f  plane  will  correspond,  in  the  case  of  a  body  wave, 
to  a  (nearly)  undispersed  signal,  with  spectral  amplitude  significantly 
above  background  in  a  frequency  range  corresponding  to  some  fraction  of 
the  total  band.  This  frequency  band  will  be  in  a  range  where  the  signal 
power  is  expected  to  be  hi ghest  relative  to  noise.  (This  means  we  will 
make  use  of  the  matched  filtering  concept,  in  the  sense  that  we  know 
roughly  what  spectral  content  we  expect  for  the  signal.  This  concept 
is  also  used  when  we  look  for  only  undispersed  signals  or  signals  of  known 
dispersion  characteristic.)  Thus,  in  the  t  -  f  plane  one  would  search 
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for  a  straight,  horizontal  line-up  of  envelope  maxima  in  a  selected  fre¬ 
quency  band,  using  the  largest  maxima  as  the  beginning  point  in  such  a 
pattern  recognition  procedure.  This,  therefore,  implicitly,  uses  a 
threshold  detection  criteria  in  a  detection  band,  since  by  starting  with  ■ 

the  largest  envelope  maxima  we  are  essentially  requiring  that  the  signal 
power  be  above  the  background  level  in  at  least  a  part  of  the  detection 
band. 

Basically,  then,  by  looking  in  this  tg  -  f  -  Ag  -  <f>g  space,  we 
will  apply  criteria  based  on  properties  of  the  expected  signal,  namely 
its  expected  dispersion  and  spectral  content,  in  order  to  recognize  a 
signal  pattern  and  to  thereby  detect  the  signal.  An  example  of  how  this 
is  accomplished  is  shown  in  Figure  2.21.  In  this  figure  we  plot  the 
times  (tg)  of  the  envelope  maxima  from  narrow  band  filter  outputs,  with  \ 

of  the  order  of  N  =  20  filters  used  so  that  the  signal  frequency  content 
is  sampled  at  about  20  points,  f  .  With  each  envelope  maxima  point  in  i 

ij 

the  plane,  there  is  also  an  associated  (spectral)  amplitude  Ag (fn ) ,  \ 

instantaneous  phase  4>g(fR)  and  an  instantaneous  frequency  (d<f>/dt)tg-  \ 

Thus,  the  t  -  f  plot  corresponds  to  a  multidimensional  display  of  ! 

9  1 

spectral  content  and  energy  arrival  time  for  a  given  segment  of  a  time  | 

A 

series.  Normally  either  1,024  or  2,048  points  are  used  for  time  seg-  I 

ments,  and  for  short  period  seismic  data  this  corresponds  to  a  50  to  100  \ 

second  segment  which  is  processed  in  each  pass.  For  on-line  continuous  | 

processing,  overlapping  time  segments  will  be  used.  j 

A  sub-band  within  the  entire  frequency  range  covered  by  the  set  j 

of  filters  is  shown  in  Figure  2.21  and  is  used  as  a  "detection  band,"  5 

that  is,  a  frequency  band  within  which  a  signal  pattern  (straight  hori-  j 

zontal  line  locus  of  envelope  maxima  in  the  case  of  an  undispersed 
body  wave)  is  sought.  This  band,  from  f£  to  f£,  is  selected  externally, 
based  on  the  expected  signal  frequency  character.  The  largest  envelope 
maxima  within  this  band  are  flagged  (denoted  by  QTj  in  the  figure)  and  \ 

used  to  compute  a  mean  "signal"  arrival  time  (i.e.,  group  time)  Tg  for  j 

the  maximum  power  arriving  in  this  frequency  range.  An  acceptance  win-  j 

dow  in  time,  for  which  maxima  can  actually  be  associated  with  an  undis-  j 

persed  signal,  can  be  constructed  using  the  relation:  t  ±  t4  ±  aAt;  ■] 

Q  0  j 

where  At  is  the  time  uncertainty  associated  with  the  envelope  maximum  f 
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Figure  2.21. 


Typical  tg-f  plane  representation  of  a  time  series  seg¬ 
ment  (0  to  36  sec)  when  the  signal-to-noise  ratio  is  low. 
Only  the  largest  envelope  maxima  (X)  and  second  largest 
maxima  (O)  are  shown.  Other,  numerous,  smaller  envelope 
maxima  are  normally  scattered  throughout  the  tg-f  plane, 
but  for  t.iis  illustration  they  have  been  omitted.  Each 
filter  "output  line"  normally  displays  about  ten  such 
peaks. 
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time  tg  for  a  filter  output  at  center  frequency  f  and  half  power  band 

width  Af.  In  particular,  AtAf  >  1/4II  is  the  theoretical  uncertainty 

relation  and  since  Af/f  =  Q’*,  then  At  _>  1/4II  (Q/f).  Further,  we  use 

t  to  denote  the  standard  deviation  in  the  time  data  used  to  compute  t„ 
o  _  9 

and  make  use  of  it  to  define  the  acceptance  time  window  about  tg  as 
well.  Hence,  taking  an  appropriately  chosen  constant  a  near  unity, 
then  two  time  window  boundary  lines  can  be  defined  to  give  the  time  win¬ 
dow;  that  is:  t  +  t*  +  (a/ 431) (Q/f)  and  t  -  t"  -  (a/4H)(Q/f) .  All 
of  the  largest  maxima  within  such  a  window  are  then  to  be  taken  as  ac¬ 
ceptable  undispersed  "signal"  group  arrivals.  Those  outside  would  be 
rejected  and  the  second  largest  maxima  could  then  be  tried,  and  so  on. 

If  the  number  of  peaks  accepted  in  the  window  is  lower  than  some 
specified  lower  limit,  then  the  tentative  "signal"  detection  would 
be  rejected  and  reclassified  as  noise." 

As  currently  implemented  for  the  SDAC  detection  experiment,  the 
algorithm  can  be  described  by  the  flow  diagram  given  in  Figure  2.22.  A 
difference  between  MARS  and  most  other  detectors  appears  right  at  the 
beginning,  in  that  quite  long  (say  100  seconds)  signal  windows  are 
typically  used. 

Fourier  transform  methods  are  used  to  derive  20  or  so  narrow  band 
envelope  functions.  A  heterodyne  operation  applied  to  the  windowed 
frequency  function  means  that  the  algorithm  effectively  calculates  one 
forward  and  two  (rather  than  20)  inverse  transforms. 

With  the  100  second  sample  of  envelope  function  available,  the 
MARS  detector  quickly  zeros  in  one  the  most  likely  candidate  "event" 
by  searching  for  the  largest  envelope  maximum  in  each  band,  and  averag¬ 
ing  their  occurrence  times.  A  generous  search  window  is  constructed 
about  this  averaged  time,  within  which  undispersed  alignments  of  envelope 
maxima  are  sought.  It  is  at  this  stage  that  a  threshold  amplitude 
value,  based  upon  the  LTA  for  each  frequency  band,  is  used  to  sort  out 
only  the  "large"  peaks.  After  this  stage,  the  actual  value  of  an  envelope 
maximum  is  never  used.  If  a  significant  number  of  frequency  bands  show 
well  aligned  envelope  peaks  (or  energy  arrivals),  an  event  is  called 
and  its  arrival  time  within  the  window  noted.  Otherwise,  the  LTA,  or 
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Figure  2.22.  MARS  detector  flow  diagram. 
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noise  threshold,  in  each  band  is  updated  and  the  next  block  of  data  pro¬ 
cessed. 

Although  a  key  part  of  the  MARS  analysis  in  the  SDAC  discrimina¬ 
tion  experiment  was  reliable  performance  of  the  detector  algorithm,  the 
dependence  of  FAR  on  program  variables  has  only  been  studied  in  the  past 
six  months.  These  noise  studies,  as  well  as  processing  of  the  SDAC 
detection  test  tape  have  yielded  quantitative  values  for  the  receiver 
operating  characteristic. 

In  the  recent  report  by  Farrell  et  al.  (1980),  it  was  concluded  that 
"The  MARS  seismic  event  detector  offers  a  significant  improvement  over  the 
current  VSC  optimally-filtered  STA/LTA  detector.  On  nearly  45  hours  of 
synthetic  data,  MARS  detected  13  percent  more  events  than  the  STA/LTA  detec¬ 
tor,  demonstrating  its  capability  of  extracting  low-level  signals  in  a  poor 
SNR  environment.  The  additional  events  detected  by  MARS  are  nearly  all 
small  ampl itude. . .events.  The  advantage  of  the  MARS  algorithm  at  low 
signal  levels  should  not  be  surprising  since  MARS  is  not  simply  a  power- 
law  detector  but  also  uses  the  signal  dispersion  and  bandwidth  as  discrim¬ 
inating  characteristics. 

The  improvement  in  detections  was  achieved  with  no  attendant  increase 
in  the  false  alarm  rate.  In  fact,  with  the  NORSAR  data,  the  MARS  FAR  was 
only  two-thirds  that  of  the  STA/LTA  detector.  The  MARS  FAR  with  the  Pine- 
dale  data  is  equal  to  that  of  the  STA/LTA  detector.  It  is  likely  that  the 
MARS  FAR  can  be  lowered  with  the  implementation  of  more  discriminating 
detection  tests,  but  the  present  level  of  MARS  detections  is  close  to  the 
theoretical  limit  predicted  for  an  ideal  matched  filter  and  is  unlikely  to 
be  greatly  increased  by  any  means. 

The  ten  percent  advantage  in  probability  of  detection  shown  by  the 
MARS  detector  for  this  particular  class  of  events  is  equivalent  to  an  im¬ 
provement  of  about  0.1  in  body  wave  magnitude.  This  is,  the  MARS  detector 
should  have  nearly  the  same  probability  of  detection  and  same  false  alarm 
rate  as  the  benchman  VSC  detector,  but  for  signals  which  are  an  average  of 
0.1  magnitude  smaller." 
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III.  CONCLUSIONS  AND  RECOMMENDATIONS 


1.  All  detectors  reviewed  in  this  report  are  based  on  a  compari¬ 
son  between  some  form  of  variance  estimate  of  the  signal  cal¬ 
culated  over  a  time  period  approximately  equal  to  the  duration 
of  an  expected  seismic  "event,"  and  the  normal  (or  long  term) 
variance. 

Possible  exceptions  to  this  general  form  are  the  phase  de¬ 
tector  described  by  Unger  (1978)  and  the  similar  part  of 
MARS.  The  former  did  not  work  because  the  short  period 
phase  fluctuations  were  so  rapid  that  they  could  not  reliably 
be  predicted;  the  latter  has  not  been  implemented.  . 

2.  Ad  hoc  detection  algorithms  such  as  Allen's  and  Stewart's 
(1977)  offer  a  hope  of  some  improvement  in  the  ROC  of 
the  detection  process.  These  methods,  after  having  made 
the  basic  detection  statistic  test  and  having  made  a  tenta¬ 
tive  identification,  are  followed  by  a  further  process  that 
is  designed  to: 

a.  reduce  the  false  alarm  rate, 

b.  improve  the  timing  capability,  and 

c.  speed  the  recovery  of  the  detection  algorithm  after  an 
event  is  encountered. 

How  much  these  algorithms  can  be  expected  to  offer  signifi¬ 
cant  improvements  in  the  ROC  depend  upon  where  on  these  curves 
they  are  operating.  If  the  SNR  is  such  that  a  high  probability 
of  detection  is  achieved  with  a  tolerable  FAR,  then  little  im¬ 
provement  can  be  expected.  If,  however,  the  SNR  is  poor  so  that 
the  probability  of  detection  is  low  for  a  given  FAR,  then  sig¬ 
nificant  improvement  may  be  achieved  by  such  two-stage 
algorithms. 

3.  Testing  of  the  detectors  described  in  this  report  was  not  car¬ 
ried  out  in  a  uniform  manner  on  the  same  or  even  similar  data 
sets.  Thus,  no  definitive  statement  about  their  relative 
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performance  may  be  made.  However,  it  appears  that  no  one 
detector  is  obviously  superior  to  all  others.  Some  have 
been  optimized  for  teleseismic  detection  using  a  moderately 
large  computing  facility.  Some  have  been  designed  for  opera¬ 
tion  in  microprocessor  controlled  field  recording  units, 
typically  used  in  the  near  field  or  at  least  at  regional 
distances.  These  applications  stress  ease  of  implementation 
in  microprocessors.  Some  algorithms  stress  detection  timing 
accuracy  as  this  may  be  another  object  of  the  algorithm. 

4.  Little  theoretical  advance  has  been  achieved  since  Frieberger1 s 
1963  work  in  which  he  "solved"  the  problem  of  detection  of  a 
Gaussian  signal  in  Gaussian  noise.  The  influence  of  this  work 
has  rightly  guided  the  design  and  implementation  of  many  of 
the  detectors  used  today.  However,  it  must  not  be  thought 
that  the  "real  problem"  has  been  solved.  The  "real  problem" 
is  the  detection  of  certain  non-Gaussian  signals  in  the  pre¬ 
sence  of  non-Gaussian  noise.  Frieberger  discovered  the  optimum 
detector  for  an  approximate  model  of  the  actual  situation, 
but  other  algorithms  may  work  significantly  better  on  "real" 
data.  Farther,  considerations  such  as  timing  ability,  re¬ 
covery  of  the  algorithm  after  an  event,  and  immunity  to  highly 
non-Gaussian  noise  such  as  line  spikes  and  data  drop  outs 
may  dictate  very  different  approaches. 

RECOMMENDATIONS 

t  Uniform  testing  of  all  viable  detectors  should  be  conducted 
with 

a.  a  realistic  synthetic  data  set, 

b.  all  detectors  coded  for  and  running  on  the  same  or 
similar  machines;  and 

c.  ROC  curves  produced  for  each  detector. 

t  Theoretical  and  experimental  research  in  phase  sensitive  de¬ 
tectors  should  be  supported  as  this  is  an  area  where  it  may 
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be  possible  to  exploit  signal  information  not  used  In 
energy  or  power  detectors. 

Research  should  be  encouraged  to  test  if  matched  filters 
for  specific  types  of  signals  (station/source  pairs)  im¬ 
proves  detector  performance. 

Further  development  of  hybrid  detection  algorithms,  that 
combine  high  probability  of  detection  with  a  high  FAR  and 
then  are  followed  by  a  post  processing  to  reduce  the  FAR 
should  be  undertaken. 

New  algorithms,  designed  specifically  for  use  with  the 
new  generation  of  three  componenet  broad  band  seismographs, 
should  be  developed  and  tested,. 

In  view  of  data  rates  from  these  new  data  sources,  dedicated 
microprocessors  to  run  the  detection  algorithms  on  a  one 
processor  per  channel  (or  seismic  station)  basis  should  be 
considered. 
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Single  Traces,"  BSSA,  68>  PP-  1521-1532. 

Abstract:  A  computer  program  has  been  developed  for  the  automatic 
detection  and  timing  of  earthquakes  on  a  single  seismic 
trace.  The  program  operates  on  line  and  is  sufficiently 
simple  that  it  is  expected  to  work  in  inexpensive  low- 
power  microprocessors  in  field  applications.  In  tests  with 
analog  tapes  of  earthquakes,  the  program  correctly  identi¬ 
fied  and  timed  to  within  0.05  sec  about  70  percent  of  the 
events  which  would  normally  be  timed  in  operation  of  a 
network.  The  program  evaluates  the  accuracy  of  its  picks, 
and  its  estimates  appear  to  be  quite  reliable.  The 
algorithm  is  working  at  present  in  a  16-bit  minicomputer 
and  appears  to  be  compatible  with  presently  available 
microprocessors. 

Farnbach,  J.  S.  (1975),  "The  Complex  Envelope  in  Seismic  Signal  Analysis," 

BSSA,  65,  pp.  951-962. 

Abstract:  Some  practical  implications  of  the  complex  envelope  repre¬ 
sentation  of  seismic  signals  are  presented.  Beginning 
with  a  look  at  an  artificially  constructed  signal  and 
proceeding  to  seismic  records,  it  is  seen  that  the  complex 
envelope  is  more  amenable  to  visual  interpretation  than 
the  real  signal  itself.  This  is  attributed  to  the  natural 
separation  of  amplitude  information  from  angle  informa¬ 
tion  afforded  by  the  complex  representation,  and  examples 
of  arrival  time  measurement  and  P-coda  correlation  suggest 
that  this  leads  to  concrete  seismological  benefits.  On 
this  basis,  it  is  suggested  that  the  complex  envelope  may 
be  a  useful  tool  in  seismic  signal  analysis. 
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Frieberger,  W.  F.  (1963),  "An  Approximate  Method  in  Signal  Detection," 
Quarterly  Appl .  Math.,  20,  pp.  373-378. 

Abstract:  A  theorem  from  the  theory  of  Tceplitz  forms  is  applied  to  the 
problem  of  estimating  the  best  test  statistic  for  the  detec¬ 
tion  of  Gaussian  signals  in  Gaussian  noise. 

Goforth,  T.  and  E.  Herrin  (1980),  "Semiannual  Technical  Report,"  Air 
Force  Office  of  Scientific  Research,  Contract  No.  F49620-76-C-0030, 
Dallas  Geophysical  Laboratory,  Southern  Methodist  University. 

Abstract:  An  automatic  seismic  signal  detection  algorithm  based  on  the 
Walsh  transform  has  been  developed.  Since  the  amplitude  of 
a  Walsh  function  is  either  +1  or  -1,  the  Walsh  transform  can 
be  accomplished  in  a  computer  with  a  series  of  shifts  and 
fixed  point  additions.  The  savings  in  computation  time  makes 
it  possible  to  compute  the  Walsh  transform  and  to  perform 
band-pass,  pre-whitening  and  adaptive  filtering  with  a  micro¬ 
computer  in  real  time  for  use  in  signal  detection. 

The  algorithm  has  been  programed  in  FORTRAN  on  a  Raytheon 
Data  Systems  500  mi ni -computer .  Tests  utilizing  seismic 
data  recorded  in  Dallas,  Albuquerque,  and  Norway  indicate 
that  the  algorithm  has  a  detection  capability  comparable  to 
a  human  analyst.  The  Walsh  detection  algorithm  runs  in 
approximately  1/10  real  time  on  the  RDS-500  mini -computer. 
Programming  of  the  detection  system  in  machine  language 
on  a  North  Star  Horizon  microprocessor-based  computer  is 
almost  complete.  Run  time  on  the  Horizon  is  estimated  to 
be  1/3  real  time. 
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Shensa,  M.  J.  (1977),  "The  Deflection  Detector,  Its  Theory  and  Evalua¬ 
tion  on  Short-Period  Seismic  Data,"  Report  TR-77-03,  Texas  Instruments, 

Alexandria,  VA. 

Abstract:  This  study  investigates  the  application  of  a  deflection  detector 
to  short-period  seismic  data.  In  general,  for  power  detectors, 
no  single  filter  will  be  optimal  for  a  large  variety  of  signals 
in  a  dynamic  noise  environment.  The  deflection  detector 
represents  an  attempt  to  adapt  to  such  a  situation  by  utilizing 
individual  FFT  frequency  cells  as  a  bank  of  filters  which  can 
accomodate  a  broad  variety  of  signals.  The  performance  of 
the  deflection  detector  is  analyzed  and  compared  to  that  of 
the  power  detector  for  several  seismic  signals.  It  is  con¬ 
cluded  that  the  deflection  detector  shows  a  distinct  advantage 
when  the  variety  of  signal  spectra  to  be  detected  is  sufficiently 
large. 

Stewart,  S.  W.  (1977),  "Real-Time  Detection  and  Location  of  Local  Seismic 
Events  in  Central  California,"  BSSA,  67,  pp.  433-452. 

Abstract:  A  computer-based  system  dedicated  full  time  to  automatic  de¬ 
tection  and  location  of  local  seismic  events  in  central 
California  has  been  developed.  The  system  monitors  108 
short-period  vertical -component  stations  from  the  U.S. 

Geological  Survey  central  California  and  Oroville  seismic 
networks.  Locations  and  matnitudes,  when  determined,  are 
printed  out  along  with  first  arrival  times,  within  2  to  5 
minutes  after  an  event  occurs.  Wave  onsets  must  be  clear 
and  impulsive  for  best  results.  For  this  reason,  regional 
events  and  teleseisms  are  usually  rejected. 

The  best  results  have  been  obtained  for  the  relatively  dense, 

16-station  Oroville  network.  For  the  month  of  October  1975, 

107  (91  percent)  of  the  118  events  timed  by  hand  were  also 
timed  and  located  by  the  real-time  system.  An  additional 
eight  events  (7  percent)  were  detected  In  real-time  but  were 
not  successfully  located.  Of  tie  107  events  for  which  both 
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on-line  and  hand-timed  locations  are  available,  92  percent 
of  the  on-line  locations  are  within  2  km  of  the  epicenters 
determined  by  hand-timing. 

During  October  1975  the  real-time  system  monitored  9i  of  the 
150  stations  of  the  central  California  network.  Of  the  260 
events  located  by  hand-timing,  225  (86  percent)  were  detected 
by  the  real-time  system.  Magnitudes  of  detected  events  range 
from  0.8  to  2.9.  Approximately  95  percent  of  the  events  of 
magnitude  U  and  greater  detected  and  located  by  hand-timing 
methods  were  also  detected  by  the  real-time  system.  Dif¬ 
ferences  between  hypocentral  locations  based  on  hand-timed 
and  computer-timed  arrivals  may  vary  from  0.1  to  5  minutes 
of  latitude  or  longitude. 

Swindell,  W.  H.  and  N.  S.  Snell  (1977),  "Station  Processor  Automatic 
Signal  Detection  System.  Phase  1:  Final  Report,  Station  Processor  Soft¬ 
ware  Development,"  Texas  Instruments  Report  No.  ALEX(Ol)-FR-77-0l, 

AFTAC  Contract  No.  F03606-76-C-0025,  Texas  Instruments,  Incorporated, 
Dallas,  Texas. 

Abstract:  This  report  summarizes  the  results  of  a  program  to  develop 
an  automatic  short-period  signal  detector  for  the  Station 
Processor  system.  Of  the  two  types  of  detectors  considered, 
the  Fisher  detector  and  the  conventional  power  detector,  the 
power  detector  was  found  to  be  superior  both  in  terms  of 
signal  response  and  false  alarm  statistics.  A  new  means  of 
setting  the  alarm  threshold  was  developed.  This  technique 
produces  a  constant  false  alarm  rate  detector  and  represents 
a  significant  improvement  over  presently  used  schemes.  A 
detection  analyzer  which  reduces  redundant  detections  from 
signal  coda  also  was  developed.  A  structure  for  the  proto¬ 
type  detection  system  was  designed  and  recommended. 
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Unger,  R.  T.  (1978),  "Automatic  Detection,  Timing  and  Preliminary  Dis¬ 
crimination  of  Seismic  Signals  with  the  Instantaneous  Amplitude,  Phase 
and  Frequency,"  Texas  Instruments  Report  No.  ALEX(01)-TR-77-04,  AFTAC 
Contract  No.  F08606-77-C-0004,  Texas  Instruments,  Incorporated,  Dallas, 
Texas. 

Abstract:  The  feasibility  Is  evaluated  of  applying  Instantaneous  ampli¬ 
tude,  phase  and  frequency  measurements  to  automatically  detect, 
time  and  Identify  seismic  events.  Detection  based  on  phase 
measurements  Is  shown  to  be  In  principle  6  dB  more  sensitive 
than  detection  based  on  amplitude  measurements.  A  phase 
detection  and  timing  algorithm,  using  a  priori  known  disper¬ 
sion  characteristics,  Is  demonstrated  to  time  the  onset  of 
simulated  teleselsmlc  long-period  surface  waves  within  30 
seconds  accuracy  In  70%  of  the  tested  cases,  for  waveforms 
down  to  0  dB  signal -to-nolse  ratio.  By  phase  measurement, 
rather  than  by  amplitude  measurement,  this  algorithm  also 
provides  a  measure  of  the  surface  wave  signal -to-nolse  ratio. 
These  results  can  be  applied  in  the  extraction  of  weak  surface 
waves . 

Phase  detection  of  teleseismic  short-period  bodywaves  was 
found  to  be  unfeasible,  due  to  the  Interference  of  early- 
arriving  secondary  signals.  Therefore,  short-period  P-wave 
detection  and  timing  are  performed  essentially  by  envelope 
peak  detection;  Instantaneous  frequency  measurements  are  also 
used  In  the  timing  process.  Tested  on  a  small  data  base, 
this  method  resulted  In  81%  to  94%  detection  at  7  to  20  false 
alarms  per  hour,  with  signal-to-noise  ratio  thresholds  of 
2  to  3  dB.  The  RMS  timing  error,  relative  to  analyst  picks, 
was  0.21  seconds,  comprising  84%  of  the  test  cases;  this 
timing  error  apparently  was  independent  of  the  signal-to- 
noise  Uio.  In  some  cases,  however,  noise  can  obscure  the 
true  signal  onset  for  the  analyst  as  well  as  for  the 
automatic  timing  algorithm.  Emergent  signals  may  cause 
timing  errors  of  several  seconds.  Measurements  of  the 
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Instantaneous  frequency  permit  analysis  of  the  delay  times 
of  secondary  signals  partially  overlapping  with  earlier 
primary  signals,  down  to  the  primary  signal  detection  level. 

Simultaneous  measurements  of  the  mean  Instantaneous  fre¬ 
quency  and  the  amount  of  Instantaneous  phase  fluctuation  over 
the  first  few  seconds  after  the  short-period  primary  signal 
onset  provided  significant  separatist  between  the  populations 
of  shallow  Eurasian  earthquakes,  Russian  presumed  nuclear 
explosions  (Including  peaceful  explosions),  and  Nevada  Test 
Site  presumed  nuclear  explosions,  even  at  signal -to-nolse 
ratios  below  0  dB. 

Vanderkulk,  W.,  F.  Rosen  and  S.  Lorenz  (1965),  "Large  Aperture  Seismic 
Array  Signal  Processing  Study,"  IBM  Final  Report,  ARPA  Contract  SD-296, 
International  Business  Machines,  Rockville,  Maryland 

Abstract:  This  report  presents  the  results  of  a  five-month  study, 

entitled  "LASA  Signal  Processing  Study"  (SD-296),  performed 
by  IBM  for  the  Advanced  Research  Projects  Agency  to: 

1.  Define  the  Large-Aperture  Seismic  Array  (LASA)  signal 
processing  requirements. 

2.  Specify  the  characteristics  of  equipments  required  to 
Implement  the  processing  requirements. 

3.  Define  an  experimental  program  to  calibrate  and  evaluate 
the  signal  processing  equipments. 

Section  1,  System  Description,  fulfills  the  requirements 
defined  under  Item  1.  Section  2,  Parametric  Analysis, 
evaluates  those  parameters  whose  physical  values  determine 
the  numerics  of  the  LASA  processing  requirements  and 
complements  Section  1. 

Section  3,  Processing  System  Configuration,  fulfills  the 
requirements  of  item  2  by  describing  the  overall  system 
in  terms  of  realizable  digital  hardware. 
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Within  Section  4,  Special  Problems,  fulfillment  of  Item  3, 
through  the  subsection  entitled  Experimental  Steering  Delay 
Determination,  Is  obtained.  In  addition,  other  topics  of 
Interest  are  Included  to  Indicate  specific  areas  of  study 
that  the  work  accomplished  has  Identified.  Within  the 
mathematical  appendixes  are  Included  the  pertinent  analytical 
studies  performed,  together  with  detailed  program  listings 
and  numerical  examples. 

Section  5,  Conclusion  and  ^commendations,  outlines  the  next 
logical  steps  In  a  program  designed  to  acquire  a  LASA  signal 
processing  capability. 
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